Fractals/fractalzoomer

Description

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)

How to run it?

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

Code repositories

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.

void createImage(double xCenter, double yCenter, double size, int image_size, int maxIterations)
{ 
   double dx = size / image_size;
   double dy = size / image_size;

   double xCenterOffset = xCenter - size * 0.5;
   double yCenterOffset = yCenter + size * 0.5;

   for(int y = 0; y < image_size; y++) {
       for(int x = 0; x < image_size; x++) {
           double result = Iterate(new Complex(xCenterOffset + dx * x, yCenterOffset - dy * y), maxIterations);
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);
           //set the color to the image at (y,x)
       }
   }
}

double Iterate(Complex value, int maxIterations)
{ 
   Complex z = value;

   for(int iterations = 0; iterations < maxIterations; iterations++)
   {
      if(triggeredCriterion(z) == true) // either escaped or converged
      {
          return outcoloring_algorithm.getResult(); //add the required arguments
      }
      
      z = function(z); // update the value based on the fractal type e.g. z = z^2 + c
   }

   return incoloring_algorithm.getResult(); //add the required arguments
}

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

double calculateFractal(Complex pixel) {

   return calculateFractalWithPeriodicity(plane.transform(rotation.rotate(pixel))); // apply plane transformation and rotation

}

double calculateFractalWithoutPeriodicity(Complex pixel) {

        int iterations = 0;

        Complex tempz = new Complex(pertur_val.getValue(init_val.getValue(pixel))); // initial value and perturbation

        Complex[] complex = new Complex[2];
        complex[0] = tempz;//z
        complex[1] = new Complex(pixel);//c

        Complex zold = new Complex(); // zold = 0, will store z(n-1)
        Complex zold2 = new Complex(); // zold2 = 0, will store z(n-2)
        Complex start = new Complex(complex[0]);
        
        Complex[] vars = createGlobalVars();

        for(; iterations < max_iterations; iterations++) {

            if(bailout_algorithm.escaped(complex[0], zold, zold2, iterations, complex[1], start, vars)) { // check if escaped
                Object[] object = {iterations, complex[0], zold, zold2, complex[1], start, vars};
                return out_color_algorithm.getResult(object);
            }
            zold2.assign(zold); // zold2 = zold
            zold.assign(complex[0]); // zold = z
            function(complex); // z = ...., for example z^2 + c

        }

        Object[] object = {complex[0], zold, zold2, complex[1], start, vars};
        return in_color_algorithm.getResult(object);
}

This is the generic code for a converging type fractal.

double calculateFractalWithoutPeriodicity(Complex pixel) {
       int iterations = 0;
       double temp = 0;

        Complex[] complex = new Complex[1];
        complex[0] = new Complex(pixel);//z

        Complex zold = new Complex(); // zold = 0, will store z(n-1)
        Complex zold2 = new Complex(); // zold2 = 0, will store z(n-2)
        Complex start = new Complex(complex[0]);
        
        Complex[] vars = createGlobalVars();

        for (; iterations < max_iterations; iterations++) {
            if((temp = complex[0].distance_squared(zold)) <= convergent_bailout) { // check if converged
                Object[] object = {iterations, complex[0], temp, zold, zold2, pixel, start, vars};
                return out_color_algorithm.getResult(object);
            }
            zold2.assign(zold); // zold2 = zold
            zold.assign(complex[0]); // zold = z
            function(complex); // z = ...., for example z - (z^3-1)/(3*z^2)
 
        }

        Object[] object = {complex[0], zold, zold2, pixel, start, vars};
        return in_color_algorithm.getResult(object);
        
}

Fractal types

Escaping Type Fractals

Fractals that check if they exceeded some predefined boundary, usually use this escape check criterion $|z_{n}|\geq bailout$ .

boolean escaped(Complex z, bouble bailout)
{ 
  if(z.norm() >= bailout)
  {
    return true;
  }

  return false;
}

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}|\leq error$ .

boolean converged(Complex z, Complex z_old, double error)
{ 
  if(z.sub(z_old).norm() <= error)
  {
    return true;
  }

  return false;
}

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

Newton

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

Halley Method

Halley

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

Schroder Method

Schroder

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

Householder Method

Householder

$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

Secant

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

Steffensen Method

Steffensen

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

Muller Method

Muller

$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^{2}f(z_{n-2})$

$B=(1+2q)f(z_{n})-(1+q)^{2}f(z_{n-1})+q^{2}f(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

Parhalley

$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

Laguerre

$z_{n+1}=z_{n}-{\frac {f(z_{n})deg}{f'(z_{n})\pm {\sqrt {(deg-1)^{2}f'(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)^{2}f'(z_{n})^{2}-deg(deg-1)f(z_{n})f''(z_{n})}}|$ and $|f'(z_{n})-{\sqrt {(deg-1)^{2}f'(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}|\geq bailout$

boolean circleCriterion(Complex z, bouble bailout)
{ 
  if(z.norm() >= bailout)
  {
    return true;
  }

  return false;
}

Square (Infinity Norm)

$max(|Re(z)|,|Im(z)|)\geq bailout$

boolean squareCriterion(Complex z, bouble bailout)
{ 
  if(Math.max(z.getAbsRe(), z.getAbsIm()) >= bailout)
  {
    return true;
  }

  return false;
}

Rhombus (One Norm)

$|Re(z)|+|Im(z)|\geq bailout$

boolean rhombusCriterion(Complex z, bouble bailout)
{ 
  if(z.getAbsRe() + z.getAbsIm() >= bailout)
  {
    return true;
  }

  return false;
}

N-Norm

$(|Re(z)|^{n}+|Im(z)|^{n})^{\frac {1}{n}}\geq bailout$

boolean nNormCriterion(Complex z, bouble bailout, double n_norm)
{ 
  if(Math.pow(Math.pow(z.getAbsRe(), n_norm) + Math.pow(z.getAbsIm(), n_norm), 1 / n_norm) >= bailout)
  {
    return true;
  }

  return false;
}

Strip

boolean stripCriterion(Complex z, bouble bailout)
{ 
  if(z.getAbsRe() >= bailout)
  {
    return true;
  }

  return false;
}

Halfplane

boolean halfplaneCriterion(Complex z, bouble bailout)
{ 
  if(z.getRe() >= bailout)
  {
    return true;
  }

  return false;
}

Field Lines

boolean fieldLinesCriterion(Complex z, Complex zold, bouble bailout)
{ 
  if(z.getRe() / zold.getRe() >= bailout&& z.getIm() / zold.getIm() >= bailout)
  {
    return true;
  }

  return false;
}

Where zold is the previous value of z.

Plane Transformation

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{aligned}e^{i\theta }z&=(\cos \theta +i\sin \theta )(x+iy)\\&=x\cos \theta +iy\cos \theta +ix\sin \theta -y\sin \theta \\&=(x\cos \theta -y\sin \theta )+i(x\sin \theta +y\cos \theta )\\&=x'+iy',\end{aligned}}$

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

Complex rotate(Complex value, double angle, Complex rotationCenter) {
        
    Complex temp = value.sub(rotationCenter); // subtract the rotation center
    Complex rotation = new Complex(Math.cos(angle), Math.sin(angle));
    return temp.times(rotation).plus(center); // multiply by the rotation and re-add the rotation center
        
}

Color

Palette

A palette is defined as a set of colors. It is usually stored in an array.

Color[] palette = new Color[N];

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.

Complex number = new Complex(2, 3);
double result = Iterate(number);
Color color = getColor(result, maxIterations);

Color getColor(double result, int maxIterations)
{
  if(result == maxIterations)
  {
     return Color.BLACK;
  }

  Color color = palette[((int)result) % N];
  return color;
}

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.

Complex number = new Complex(2, 3);
double result = Iterate(number);
Color color = getColor(result, maxIterations);

Color getColor(double result, int maxIterations)
{
  if(result == maxIterations)
  {
     return Color.BLACK;
  }

  Color color1 = palette[((int)result) % N];
  Color color1 = palette[((int)result + 1) % N];

  double fractional = result - (int)result;

  Color finalColor = new Color((int)(color1.getRed() + (color2.getRed() - color1.getRed()) * fractional + 0.5),
                             (int)(color1.getGreen() + (color2.getGreen() - color1.getGreen()) * fractional + 0.5),
                             (int)(color1.getBlue() + (color2.getBlue() - color1.getBlue()) * fractional + 0.5));
  return finalColor;
}

Algorithms

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.

Maximum Iterations

double maximumIterations(int maxIterations)
{ 
   return maxIterations; //max iterations
}

norm(z)

double normZ(Complex z, int maxIterations)
{ 
   return z.norm_squared() * (maxIterations / 3.0);
}

Decomposition Like

double decompositionLike(Complex z)
{ 
   return Math.abs((z.arg() / 2 * Math.PI + 0.75) * 59 * Math.PI);
}

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

Re / Im

double ReDivideIm(Complex z)
{ 
   return Math.abs(z.getRe() / z.getIm()) * 8;
}

cos(norm(z))

double CosNormZ(Complex z)
{ 
   if((int)(z.norm_squared() * 10) % 2 == 1)
   {
      return Math.abs(Math.cos(z.getRe() * z.getIm() * z.getAbsRe() * z.getAbsIm()) * 400 + 50;
   }

   return Math.abs(Math.sin(z.getRe() * z.getIm() * z.getAbsRe() * z.getAbsIm())) * 400;
}

norm(z) * cos(Re^2)

double NormZCosReSquared(Complex z)
{ 
   return z.norm_squared() * Math.abs(Math.cos(z.getRe() * z.getRe())) * 400;
}

sin(Re^2 - Im^2)

double SinReSquaredMinusImSquared(Complex z)
{ 
   return Math.abs(Math.sin(z.getRe() * z.getRe() - z.getIm() * z.getIm())) * 400;
}

atan(Re * Im * |Re| * |Im|)

double AtanReTimesImTimesAbsReTimesAbsIm(Complex z)
{ 
   return Math.abs(Math.atan(z.getRe() * z.getIm() * z.getAbsRe() * z.getAbsIm())) * 400;
}

Squares

double Squares(Complex z)
{ 
   if(((Math.abs((int)(z.getRe() * 40)) % 2) ^ (Math.abs((int)(z.getIm() * 40)) % 2)) == 1)
   {
      return Math.abs((Math.atan2(z.getIm(), z.getRe()) / (Math.PI * 2)  + 0.75) * Math.PI * 59);
   }
   
   return Math.abs((Math.atan2(z.getRe(), z.getIm() / (Math.PI * 2)  + 0.75) * Math.PI * 59);
}

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

Squares 2

double Squares2(Complex z)
{ 
   double x = z.getRe() * 16;
   double y = z.getIm() * 16;
        
   double dx = Math.abs(x - Math.floor(x));
   double dy = Math.abs(y - Math.floor(y));
   
   if((dx < 0.5 && dy < 0.5) || (dx > 0.5 && dy > 0.5))
   {
       return 50; // a palette offset
   }

   return 0;
}

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.

double escapeTime(int n)
{ 
  return n; //iterations
}

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

Method 1

 $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 2

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

Increasing the bailout value will produce smoother images.

Calculating Continuous Iteration Count for Converging Type Fractals

Method 1

 $error=1e-9$

 $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})}}$

Method 2

 $error=1e-9$

 $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$ .

Method 1

 $error=1e-9$

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

Method 2

 $error=1e-9$

 $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.

double BinaryDecomposition(int n, Complex z)
{ 
   if(z.getIm() < 0)
   {
     return n + 50; //You can add A or B or C or D to include Continuous Iteration Count 
   }
   
   return n; //You can add A or B or C or D to include Continuous Iteration Count 
}

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.

double BinaryDecomposition2(int n, Complex z)
{ 
   if(z.getRe() < 0)
   {
     return n + 50; //You can add A or B or C or D to include Continuous Iteration Count 
   }
   
   return n; //You can add A or B or C or D to include Continuous Iteration Count 
}

Escape Time + Re

double EscapeTimePlusRe(int n, Complex z)
{ 
  return n + z.getRe();
}

Escape Time + Im

double EscapeTimePlusIm(int n, Complex z)
{ 
  return n + z.getIm();
}

Escape Time + Re/Im

double EscapeTimePlusReDivideIm(int n, Complex z)
{ 
  return n + z.getRe() / z.getIm();
}

Escape Time + Re + Im + Re/Im

double EscapeTimePlusRePlusImPlusReDivideIm(int n, Complex z)
{ 
  return n + z.getRe() + z.getIm() + z.getRe() / z.getIm();
}

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.

double Biomorph(int n, Complex z, double bailout)
{ 
   if(z.getRe() > -bailout && z.getRe() < bailout
   || z.getIm() > -bailout && z.getIm() < bailout)
   {
      return n; //You can add A or B to include Continuous Iteration Count
   }
   
   return n + 50; //You can add A or B to include Continuous Iteration Count 
}

Color Decomposition

double ColorDecomposition(int n, Complex z)
{ 
  return Math.abs((z.arg() / 2 * Math.PI + 0.75) * 59 * Math.PI);
}

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.

double ColorDecomposition(int n, Complex z)
{ 
  double re = Math.floor(1000 * (z.getRe() + 0.5) / 1000; // rounding
  double im = Math.floor(1000 * (z.getRe() + 0.5) / 1000; // rounding
  
  return Math.abs((Math.atan2(im, re) / Math.PI * 2  + 0.75) * Math.PI * 59  + (re * re + im * im) * 2.5);
}

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

Escape Time + Color Decomposition

double EscapeTimeColorDecomposition(int n, Complex z)
{ 
  return n + Math.abs((z.arg() / 2 * Math.PI + 0.75) * 59 * Math.PI);
}

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.

double EscapeTimeColorDecomposition(int n, Complex z)
{ 
  double re = Math.floor(1000 * (z.getRe() + 0.5) / 1000; // rounding
  double im = Math.floor(1000 * (z.getRe() + 0.5) / 1000; // rounding
  
  return n + Math.abs((Math.atan2(im, re) / Math.PI * 2  + 0.75) * Math.PI * 59  + (re * re + im * im) * 2.5); //You can add C or D to include Continuous Iteration Count 
}

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

Escape Time + Gaussian Integer

double EscapeTimeGaussianInteger(int n, Complex z)
{ 
  return n + z.sub(z.gaussian_integer()).norm_squared() * 90;
}

Where gaussian integer can be obtained with this function:

Complex gaussian_integer()
{ 
  return new Complex((int)(re < 0 ? re - 0.5 : re + 0.5), (int)(im < 0 ? im - 0.5 : im + 0.5));
}

Escape Time + Gaussian Integer 2

double EscapeTimeGaussianInteger2(int n, Complex z)
{ 
  Complex temp = z.sub(z.gaussian_integer());
  return n + Math.abs(Math.atan(temp.getIm() / temp.getRe())) * 5;
}

Escape Time + Gaussian Integer 3

double EscapeTimeGaussianInteger3(int n, Complex z)
{ 
  Complex temp = z.sub(z.gaussian_integer());
  return Math.abs(n + temp.getRe());
}

Escape Time + Gaussian Integer 4

double EscapeTimeGaussianInteger4(int n, Complex z)
{ 
  Complex temp = z.sub(z.gaussian_integer());
  return Math.abs(n +  temp.getRe() + temp.getIm());
}

Escape Time + Gaussian Integer 5

double EscapeTimeGaussianInteger5(int n, Complex z)
{ 
  Complex temp = z.sub(z.gaussian_integer());
  return Math.abs(n +  temp.getRe() + temp.getIm() + temp.getRe() / temp.getIm());
}

Escape Time + Algorithm

double EscapeTimePlusAlgorithm(int n, Complex z, Complex z_old)
{ 
  Complex temp = z.sub(z_old);
  return n + Math.abs(Math.atan(temp.getIm() / temp.getRe())) * 4;
}

Escape Time + Algorithm 2

double EscapeTimePlusAlgorithm2(int n, Complex z)
{ 
  Complex temp = z.sub(z.sin());
  return n + Math.abs(Math.atan(temp.getIm() / temp.getRe())) * 8;
}

Distance Estimator

Escape Time + Escape Radius

double EscapeTimeEscapeRadius(int n, Complex z, double bailout)
{ 
  double zabs = Math.log(z.norm_squared()) / Math.log(bailout * bailout) - 1.0;
  double zarg = (z.arg() / (2 * Math.PI) + 1.0) % 1.0;
  return n + zabs + zarg;
}

Escape Time + Grid

double EscapeTimeGrid(int n, Complex z, double bailout)
{ 
  double zabs = Math.log(z.norm_squared()) / Math.log(bailout * bailout) - 1.0;
  double zarg = (z.arg() / (2 * Math.PI) + 1.0) % 1.0;
  
  if(0.05 < zabs && zabs < 0.95 && 0.05 < zarg && zarg < 0.95)
  {
    return n; //You can add A or B to include Continuous Iteration Count
  }

  return n + 50; //You can add A or B to include Continuous Iteration Count
}

For smoother grid lines the following version can be used.

double EscapeTimeGrid(int n, Complex z, double bailout)
{ 
  double zabs = Math.log(z.norm_squared()) / Math.log(bailout * bailout) - 1.0;
  double zarg = (z.arg() / (2 * Math.PI) + 1.0) % 1.0;

  double k = Math.pow(0.5, 0.5 - zabs);
  double grid_weight = 0.05;
  
  if(grid_weight < zabs && zabs < (1.0 - grid_weight) && (grid_weight * k) < zarg && zarg < (1.0 - grid_weight * k))
  {
    return n; //You can add A or B to include Continuous Iteration Count
  }

  return n + 50; //You can add A or B to include Continuous Iteration Count
}

Banded

double Banded(int n, Complex z)
{ 
  return n + Math.abs((Math.log(Math.log(z.norm_squared())) / Math.log(2)) * 2.4);
}

Escape Time + Field Lines

double EscapeTimeFieldLines(int n, Complex z, double bailout)
{ 
  double lineWidth = 0.008;
  double fx = (z.arg() / (2 * Math.PI);  // angle within cell
  double fy = Math.log(z.norm_squared()) / Math.log(bailout * bailout);
  double fz = Math.pow(0.5, -fy);

  if(Math.abs(fx) > lineWidth * fz)
  {
    return n; //You can add A or B to include Continuous Iteration Count
  }
                
  return n + 50; //You can add A or B to include Continuous Iteration Count
}

Escape Time + Field Lines 2

double EscapeTimeFieldLines2(int n, Complex z, double bailout)
{ 
  double lineWidth = 0.07;
  double fx = (z.arg() / 2) * Math.PI;
  double fy = Math.log(z.norm_squared()) / Math.log(bailout * bailout);
  double fz = Math.pow(0.5, -fy);

  if(Math.abs(fx) < (1.0 - lineWidth) * fz && lineWidth * fz < Math.abs(fx))
  {
    return n; //You can add A or B to include Continuous Iteration Count
  }
                
  return n + 50; //You can add A or B to include Continuous Iteration Count
}

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."^[1]

Polar Projection

Without polar projection

With polar projection

Polar Projection or Exponential Mapping is a technique that maps the complex plane into radius and angle.

void createPolarProjectionImage(double xCenter, double yCenter, double size, int image_size, int maxIterations)
{ 
    double center = Math.log(size);
    double dy = (2 * Math.PI) / image_size;
    double dx = dy;
    double start = center - dx * image_size * 0.5;

   for(int y = 0; y < image_size; y++) {
       for(int x = 0; x < image_size; x++) {         
           double sf = Math.sin(y * dy);
           double cf = Math.cos(y * dy);
           double r = Math.exp(x * dx + start);

           double result = Iterate(new Complex(xCenter + r * cf, yCenter + r * sf), maxIterations);
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);
           //set the color to the image at (y,x)
       }
   }
}

Optimizations

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...

void createImage(int gridI, int gridJ, int gridSize, double xCenter, double yCenter, double size, int image_size, int maxIterations)
{ 
   double dx = size / image_size;
   double dy = size / image_size;

   double xCenterOffset = xCenter - size * 0.5;
   double yCenterOffset = yCenter + size * 0.5;

   int fromY = gridI * image_size / gridSize;
   int toY = (gridI + 1) * image_size / gridSize;

   int fromX = gridJ * image_size / gridSize;
   int toX = (gridJ + 1) * image_size / gridSize;

   for(int y = fromY; y < toY; y++) {
       for(int x = fromX; x < toX; x++) {
           double result = Iterate(new Complex(xCenterOffset + dx * x, yCenterOffset - dy * y), maxIterations);
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);
           //set the color to the image at (y,x)
       }
   }
}

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

void createPolarProjectionImage(int gridI, int gridJ, int gridSize, double xCenter, double yCenter, double size, int image_size, int maxIterations)
{ 
    double center = Math.log(size);
    double dy = (2 * Math.PI) / image_size;
    double dx = dy;
    double start = center - dx * image_size * 0.5;

    int fromY = gridI * image_size / gridSize;
    int toY = (gridI + 1) * image_size / gridSize;

    int fromX = gridJ * image_size / gridSize;
    int toX = (gridJ + 1) * image_size / gridSize;

   for(int y = fromY; y < toY; y++) {
       for(int x = fromX; x < toX; x++) {         
           double sf = Math.sin(y * dy);
           double cf = Math.cos(y * dy);
           double r = Math.exp(x * dx + start);

           double result = Iterate(new Complex(xCenter + r * cf, yCenter + r * sf), maxIterations);
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);
           //set the color to the image at (y,x)
       }
   }
}

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

Boundary Tracing, multithreaded

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

Divide And Conquer

Divide And Conquer, iteration visualization

Divide And Conquer, multithreaded

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.

Complex period = new Complex(); // period = 0
int check = 3;
int check_counter = 0;

int update = 10;
int update_counter = 0;

boolean periodicityCheck(Complex z) {

        //Check for period
        if(z.distance_squared(period) < 1e-13) { // |z-period|^2
            return true;
        }

        //Update history
        if(check == check_counter) {
            check_counter = 0;

            //Double the value of check
            if(update == update_counter) {
                update_counter = 0;
                check <<= 1;
            }
            update_counter++;

            period.assign(z); // period = z
        } //End of update history

        check_counter++;

        return false;
}

Super Sampling

Supersampling^[2] 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.

Super Sampling Method

This method also includes the Threaded optimization.

void createImageSuperSampling(int gridI, int gridJ, int gridSize, double xCenter, double yCenter, double size, int image_size, int maxIterations, int samples)
{ 
   double dx = size / image_size;
   double dy = size / image_size;

   double xCenterOffset = xCenter - size * 0.5;
   double yCenterOffset = yCenter + size * 0.5;

   int fromY = gridI * image_size / gridSize;
   int toY = (gridI + 1) * image_size / gridSize;

   int fromX = gridJ * image_size / gridSize;
   int toX = (gridJ + 1) * image_size / gridSize;

   double a = size * 0.25;

   double totalSamples = samples + 1;

   double x[] = {-a, a, a, -a,
            -a, a, 0, 0,
            -2*a, -2*a, -2*a, 0, 0, 2*a, 2*a, 2*a,
            -2*a, -2*a, -a, -a, a, a, 2*a, 2*a};
   double y[] = {-a, -a, a, a,
            0, 0, -a, a,
            -2*a, 0, 2*a, -2*a, 2*a, -2*a, 0, 2*a,
            -a, a, -2*a, 2*a, -2*a, 2*a, -a, a};

   for(int y = fromY; y < toY; y++) {
       for(int x = fromX; x < toX; x++) {
           double x0 = xCenterOffset + dx * x;
           double y0 = yCenterOffset - dy * y;
           double result = Iterate(new Complex(x0, y0), maxIterations); //calculate for the center value
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);
           
           int red = color.getRed();
           int green = color.getGreen();
           int blue = color.getBlue();

           for(int i = 0; i < samples; i++)
           {
               result = Iterate(new Complex(x0 + x[i], y0 + y[i]), maxIterations); //calculate for the extra samples around center
               color = getColor(result, maxIterations);

               /* Sum all the samples */
               red += color.getRed();
               green += color.getGreen();
               blue += color.getBlue();
           }

           Color finalColor = new Color((int)(red / totalSamples + 0.5),
                                        (int)(green / totalSamples + 0.5),
                                        (int)(blue / totalSamples + 0.5)); //average the color

           //set the final color to the image at (y,x)
       }
   }
}

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

void createPolarProjectionImageSuperSampling(int gridI, int gridJ, int gridSize, double xCenter, double yCenter, double size, int image_size, int maxIterations, int samples)
{ 
    double center = Math.log(size);
    double dy = (2 * Math.PI) / image_size;
    double dx = dy;
    double start = center - dx * image_size * 0.5;

    int fromY = gridI * image_size / gridSize;
    int toY = (gridI + 1) * image_size / gridSize;

    int fromX = gridJ * image_size / gridSize;
    int toX = (gridJ + 1) * image_size / gridSize;

    double a = dy * 0.25;

    double totalSamples = samples + 1;

    double x[] = {-a, a, a, -a,
            -a, a, 0, 0,
            -2*a, -2*a, -2*a, 0, 0, 2*a, 2*a, 2*a,
            -2*a, -2*a, -a, -a, a, a, 2*a, 2*a};
    double y[] = {-a, -a, a, a,
            0, 0, -a, a,
            -2*a, 0, 2*a, -2*a, 2*a, -2*a, 0, 2*a,
            -a, a, -2*a, 2*a, -2*a, 2*a, -a, a};

    for(int y = fromY; y < toY; y++) {
       for(int x = fromX; x < toX; x++) {         
           double sf = Math.sin(y * dy);
           double cf = Math.cos(y * dy);
           double r = Math.exp(x * dx + start);

           double result = Iterate(new Complex(xCenter + r * cf, yCenter + r * sf), maxIterations);
           imageIterations[y][x] = result; //save the iteration data for later use
           Color color = getColor(result, maxIterations);

           int red = color.getRed();
           int green = color.getGreen();
           int blue = color.getBlue();

           for(int i = 0; i < samples; i++)
           {
               sf = Math.sin(y * dy + y[i]);
               cf = Math.cos(y * dy + y[i]);
               r = Math.exp(x * dx + start + x[i]);

               result = Iterate(new Complex(xCenter + r * cf, yCenter + r * sf), maxIterations); //calculate for the extra samples around center
               color = getColor(result, maxIterations);

               /* Sum all the samples */
               red += color.getRed();
               green += color.getGreen();
               blue += color.getBlue();
           }

           Color finalColor = new Color((int)(red / totalSamples + 0.5),
                                        (int)(green / totalSamples + 0.5),
                                        (int)(blue / totalSamples + 0.5)); //average the color
           //set the color to the image at (y,x)
       }
    }
}

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

jar file = program ( executable file)
*.fzs file Fractal Zoomer Settings ( binary file)

References

This article is issued from Wikibooks. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[1] ractalforums: entropy-mandelbrot-coloring

[2] wikipedia: Supersampling