Fractals/Iterations in the complex plane/Mandelbrot set/mandelbrot
Speed improvements - optimisation
Symmetry Edit
The Mandelbrot set is symmetric with respect to the x-axis in the plane :
colour(x,y) = colour(x,-y)
its intersection with the x-axis ( real slice of Mandelbrot set ) is an interval :
<-2 ; 1/4>
It can be used to speed up computations ( up to 2-times )[1]
Bailout test Edit
Instead of checking if magnitude ( radius = abs(z) ) is greater then escape radius ( ER):
compute square of ER:
ER2 = ER*ER
once and check :
It gives the same result and is faster.
interior detection Edit
Interior detection method[2]
- time with detection versus without detection is : 0.383s versus 8.692s so it is 23 times faster !!!!
// gives last iterate = escape time
// output 0< i < iMax
int iterate(double complex C , int iMax)
{
int i=0;
double complex Z= C; // initial value for iteration Z0
complex double D = 1.0; // derivative with respect to z
for(i=0;i<iMax;i++)
{ if(cabs(Z)>EscapeRadius) break; // exterior
if(cabs(D)<eps) break; // interior
D = 2.0*D*Z;
Z=Z*Z+C; // complex quadratic polynomial
}
return i;
}
Period detection Edit
Period of Complex_quadratic_map - main article
"When rendering a Mandelbrot or Julia set, the most time-consuming parts of the image are the “black areas”. In these areas, the iterations never escape to infinity, so every pixel must be iterated to the maximum limit. Therefore, much time can be saved by using an algorithm to detect these areas in advance, so that they can be immediately coloured black, rather than rendering them in the normal way, pixel by pixel. FractalNet uses a recursive algorithm to split the image up into blocks, and tests each block to see whether it lies inside a “black area”. In this way, large areas of the image can be quickly coloured black, often saving a lot of rendering time. ... (some) blocks were detected as “black areas” and coloured black immediately, without having to be rendered. (Other) blocks were rendered in the normal way, pixel by pixel." Michael Hogg [3]
// cpp code by Geek3
// http://commons.wikimedia.org/wiki/File:Mandelbrot_set_rainbow_colors.png
bool outcircle(double center_x, double center_y, double r, double x, double y)
{ // checks if (x,y) is outside the circle around (center_x,center_y) with radius r
x -= center_x;
y -= center_y;
if (x * x + y * y > r * r)
return(true);
return(false);
// skip values we know they are inside
if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
&& outcircle(-1.0, 0.0, 0.25, x0, y0)
&& outcircle(-0.125, 0.744, 0.092, x0, y0)
&& outcircle(-1.308, 0.0, 0.058, x0, y0)
&& outcircle(0.0, 0.25, 0.35, x0, y0))
{
// code for iteration
}
Cardioid and period-2 checking Edit
One way to improve calculations is to find out beforehand whether the given point lies within the cardioid or in the period-2 bulb. Before passing the complex value through the escape time algorithm, first check if:
to determine if the point lies within the period-2 bulb and
to determine if the point lies inside the main cardioid.
another description [4]
// http://www.fractalforums.com/new-theories-and-research/quick-(non-iterative)-rejection-filter-for-mandelbrotbuddhabrot-with-benchmark/
public static void quickRejectionFilter(BlockingCollection<Complex> input, BlockingCollection<Complex> output)
{
foreach(var item in input.GetConsumingEnumerable())
{
if ((Complex.Abs(1.0 - Complex.Sqrt(Complex.One - (4 * item))) < 1.0)) continue;
if (((Complex.Abs(item - new Complex(-1, 0))) < 0.25)) continue;
if ((((item.Real + 1.309) * (item.Real + 1.309)) + item.Imaginary * item.Imaginary) < 0.00345) continue;
if ((((item.Real + 0.125) * (item.Real + 0.125)) + (item.Imaginary - 0.744) * (item.Imaginary - 0.744)) < 0.0088) continue;
if ((((item.Real + 0.125) * (item.Real + 0.125)) + (item.Imaginary + 0.744) * (item.Imaginary + 0.744)) < 0.0088) continue;
//We tried every known quick filter and didn't reject the item, adding it to next queue.
output.Add(item);
}
}
See also
Periodicity checking Edit
Most points inside the Mandelbrot set oscillate within a fixed orbit. There could be anything from ten to tens of thousands of points in between, but if an orbit ever reaches a point where it has been before then it will continually follow this path, will never tend towards infinity and hence is in the Mandelbrot set.
This Mandelbrot processor includes:
- periodicity checking
- period-2 bulb and cardioid checking
for a great speed up during deep zoom animations with a high maximum iteration value.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Mandelbrot2
{
public struct MandelbrotData
{
private double px;
private double py;
private double zx;
private double zy;
private long iteration;
private bool inSet;
private HowFound found;
public MandelbrotData(double px, double py)
{
this.px = px;
this.py = py;
this.zx = 0.0;
this.zy = 0.0;
this.iteration = 0L;
this.inSet = false;
this.found = HowFound.Not;
}
public MandelbrotData(double px, double py,
double zx, double zy,
long iteration,
bool inSet,
HowFound found)
{
this.px = px;
this.py = py;
this.zx = zx;
this.zy = zy;
this.iteration = iteration;
this.inSet = inSet;
this.found = found;
}
public double PX
{
get { return this.px; }
}
public double PY
{
get { return this.py; }
}
public double ZX
{
get { return this.zx; }
}
public double ZY
{
get { return this.zy; }
}
public long Iteration
{
get { return this.iteration; }
}
public bool InSet
{
get { return this.inSet; }
}
public HowFound Found
{
get { return this.found; }
}
}
public enum HowFound { Max, Period, Cardioid, Bulb, Not }
class MandelbrotProcess
{
private long maxIteration;
private double bailout;
public MandelbrotProcess(long maxIteration, double bailout)
{
this.maxIteration = maxIteration;
this.bailout = bailout;
}
public MandelbrotData Process(MandelbrotData data)
{
double zx;
double zy;
double xx;
double yy;
double px = data.PX;
double py = data.PY;
yy = py * py;
#region Cardioid check
//Cardioid
double temp = px - 0.25;
double q = temp * temp + yy;
double a = q * (q + temp);
double b = 0.25 * yy;
if (a < b)
return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Cardioid);
#endregion
#region Period-2 bulb check
//Period-2 bulb
temp = px + 1.0;
temp = temp * temp + yy;
if (temp < 0.0625)
return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Bulb);
#endregion
zx = px;
zy = py;
int check = 3;
int checkCounter = 0;
int update = 10;
int updateCounter = 0;
double hx = 0.0;
double hy = 0.0;
for (long i = 1; i <= this.maxIteration; i++)
{
//Calculate squares
xx = zx * zx;
yy = zy * zy;
#region Bailout check
//Check bailout
if (xx + yy > this.bailout)
return new MandelbrotData(px, py, zx, zy, i, false, HowFound.Not);
#endregion
//Iterate
zy = 2.0 * zx * zy + py;
zx = xx - yy + px;
#region Periodicity check
//Check for period
double xDiff = Math.Abs(zx - hx);
if (xDiff < this.ZERO)
{
double yDiff = Math.Abs(zy - hy);
if (yDiff < this.ZERO)
return new MandelbrotData(px, py, zx, zy, i, true, HowFound.Period);
} //End of on zero tests
//Update history
if (check == checkCounter)
{
checkCounter = 0;
//Double the value of check
if (update == updateCounter)
{
updateCounter = 0;
check *= 2;
}
updateCounter++;
hx = zx;
hy = zy;
} //End of update history
checkCounter++;
#endregion
} //End of iterate for
#region Max iteration
return new MandelbrotData(px, py, zx, zy, this.maxIteration, true, HowFound.Max);
#endregion
}
private double ZERO = 1e-17;
}
}
References Edit
- T Myers: Introduction to perturbation ...
- interpretation-of-topology-optimization
- Derivative Extrapolation From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2023.