View source for Plotting algorithms for the Mandelbrot set

{{short description|Algorithms and methods of plotting the Mandelbrot set on a computing device}}
{{tone|reads like a textbook|date=July 2021}}
{{Use dmy dates|date=February 2020}}
[[File:Fractal-zoom-1-03-Mandelbrot Buzzsaw.png|right|thumbnail|Still image of [//upload.wikimedia.org/wikipedia/commons/0/07/Fractal-zoom-1-03-Mandelbrot_Buzzsaw.ogv a movie of increasing magnification] on 0.001643721971153 − 0.822467633298876''i'']]
[[File:Mandelbrot sequence new still.png|right|thumbnail|Still image of [//upload.wikimedia.org/wikipedia/commons/5/51/Mandelbrot_sequence_new.webm an animation of increasing magnification]]]
There are many programs and [[algorithm]]s used to <!--generate-->plot the [[Mandelbrot set]] and other [[fractal]]s, some of which are described in [[fractal-generating software]]. These programs use a variety of algorithms to determine the color of individual [[pixel]]s efficiently. <!-- There have been many algorithms developed to efficiently plot the Mandelbrot set via a computing device. Some of the more prominent ones are discussed below. -->

==Escape time algorithm==

The simplest algorithm for generating a representation of the Mandelbrot set is known as the "escape time" algorithm. A repeating calculation is performed for each ''x'', ''y'' point in the plot area and based on the behavior of that calculation, a color is chosen for that pixel.

===Unoptimized naïve escape time algorithm===
<!--[[Spatial anti-aliasing]] links directly here.-->

In both the unoptimized and optimized escape time algorithms, the ''x'' and ''y'' locations of each point are used as starting values in a repeating, or iterating calculation (described in detail below). The result of each iteration is used as the starting values for the next. The values are checked during each iteration to see whether they have reached a critical "escape" condition, or "bailout". If that condition is reached, the calculation is stopped, the pixel is drawn, and the next ''x'', ''y'' point is examined. For some starting values, escape occurs quickly, after only a small number of iterations. For starting values very close to but not in the set, it may take hundreds or thousands of iterations to escape. For values within the Mandelbrot set, escape will never occur. The programmer or user must choose how many iterations–or how much "depth"–they wish to examine. The higher the maximal number of iterations, the more detail and subtlety emerge in the final image, but the longer time it will take to calculate the fractal image.

Escape conditions can be simple or complex. Because no [[complex number]] with a real or imaginary part greater than 2 can be part of the set, a common bailout is to escape when either coefficient exceeds 2. A more computationally complex method that detects escapes sooner, is to compute distance from the origin using the [[Pythagorean theorem]], i.e., to determine the [[absolute value]], or ''modulus'', of the complex number. If this value exceeds 2, or equivalently, when the sum of the squares of the real and imaginary parts exceed 4, the point has reached escape. More computationally intensive rendering variations include the [[Buddhabrot]] method, which finds escaping points and plots their iterated coordinates.

The color of each point represents how quickly the values reached the escape point. Often black is used to show values that fail to escape before the iteration limit, and gradually brighter colors are used for points that escape. This gives a visual representation of how many cycles were required before reaching the escape condition.

To render such an image, the region of the complex plane we are considering is subdivided into a certain number of [[pixel]]s. To color any such pixel, let <math>c</math> be the midpoint of that pixel. We now iterate the critical point 0 under <math>P_c</math>, checking at each step whether the orbit point has modulus larger than 2. When this is the case, we know that <math>c</math> does not belong to the Mandelbrot set, and we color our pixel according to the number of iterations used to find out. Otherwise, we keep iterating up to a fixed number of steps, after which we decide that our parameter is "probably" in the Mandelbrot set, or at least very close to it, and color the pixel black.

In [[pseudocode]], this algorithm would look as follows. The algorithm does not use complex numbers and manually simulates complex-number operations using two real numbers, for those who do not have a [[complex data type]]. The program may be simplified if the programming language includes complex-data-type operations.

<!-- NOTE that xtemp is necessary, otherwise y would be calculated with the new x, which would be wrong. Also note that one must plot (''x''<sub>0</sub>,&nbsp;''y''<sub>0</sub>), not (''x'',''y''). -->
 '''for each''' pixel (Px, Py) on the screen '''do'''
     x0 := scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.00, 0.47))
     y0 := scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1.12, 1.12))
     x := 0.0
     y := 0.0
     iteration := 0
     max_iteration := 1000
     '''while''' (x*x + y*y ≤ 2*2 AND iteration < max_iteration) '''do'''
         xtemp := x*x - y*y + x0
         y := 2*x*y + y0
         x := xtemp
         iteration := iteration + 1
  
     color := palette[iteration]
     plot(Px, Py, color)

Here, relating the pseudocode to <math>c</math>, <math>z</math> and <math>P_c</math>:
* <math> z = x + iy\ </math>
* <math> z^2 = x^2 +2ixy</math> - <math>y^2 \ </math>
* <math> c = x_0 + i y_0\ </math>
and so, as can be seen in the pseudocode in the computation of ''x'' and ''y'':
* <math>x = \mathop{\mathrm{Re}}(z^2+c) = x^2-y^2 + x_0</math> and <math>y = \mathop{\mathrm{Im}}(z^2+c) = 2xy + y_0.\ </math>

To get colorful images of the set, the assignment of a color to each value of the number of executed iterations can be made using one of a variety of functions (linear, exponential, etc.). One practical way, without slowing down calculations, is to use the number of executed iterations as an entry to a [[Palette (computing)|palette]] initialized at startup. If the color table has, for instance, 500 entries, then the color selection is ''n''&nbsp;mod&nbsp;500, where ''n'' is the number of iterations.

===Optimized escape time algorithms===
The code in the previous section uses an unoptimized inner while loop for clarity. In the unoptimized version, one must perform five multiplications per iteration. To reduce the number of multiplications the following code for the inner while loop may be used instead:

 x2:= 0
 y2:= 0
 w:= 0
 
 '''while''' (x2 + y2 ≤ 4 '''and''' iteration < max_iteration) '''do'''
     x:= x2 - y2 + x0
     y:= w - x2 - y2 + y0
     x2:= x * x
     y2:= y * y
     w:= (x + y) * (x + y)
     iteration:= iteration + 1

The above code works via some algebraic simplification of the complex multiplication:

<math>\begin{align}
(iy + x)^2 &= -y^2 + 2iyx + x^2 \\
 &= x^2 - y^2 + 2iyx
\end{align} </math>

Using the above identity, the number of multiplications can be reduced to three instead of five.

The above inner while loop can be further optimized by expanding ''w'' to

<!-- <math>w = xx + 2yx + yy</math> -->
<math>w = x^{2} + 2xy + y^{2}</math>

Substituting ''w'' into <math>y = w - x^{2} - y^{2} + y_0</math> yields
<math>y = 2xy + y_0</math>
and hence calculating ''w'' is no longer needed.

The further optimized pseudocode for the above is:

 x:= 0
 y:= 0
 x2:= 0
 y2:= 0
 
 '''while''' (x2 + y2 ≤ 4 '''and''' iteration < max_iteration) '''do'''
     y:= 2 * x * y + y0
     x:= x2 - y2 + x0
     x2:= x * x
     y2:= y * y
     iteration:= iteration + 1

Note that in the above pseudocode, <math>2xy</math> seems to increase the number of multiplications by 1, but since 2 is the multiplier the code can be optimized via <math>(x + x)y</math>.

=== Derivative Bailout or "derbail" ===
[[File:Derbail method render.png|thumb|379x379px|An example of the fine detail possible with the usage of derbail, rendered with 1024 samples]]
It is common to check the magnitude of z after every iteration, but there is another method we can use that can converge faster and reveal structure within [[julia set]]s.

Instead of checking if the [[Magnitude of a complex number|magnitude of z]] after every iteration is larger than a given value, we can instead check if the sum of each [[derivative]] of z up to the current iteration step is larger than a given bailout value{{Citation needed|date=September 2024}}:

<math>z_n^\prime := (2 * z_{n-1}^\prime * z_{n-1}) + 1 </math>

The size of the dbail value can enhance the detail in the structures revealed within the dbail method with very large values.

It is possible to find derivatives automatically by leveraging [[Automatic differentiation]] and computing the iterations using [[Dual number]]s{{Citation needed|date=September 2024}}. 
[[File:Example of derbail precision issues.png|thumb|253x253px|Hole caused by precision issues]]
Rendering fractals with the derbail technique can often require a large number of samples per pixel, as there can be [[Floating-point arithmetic|precision]] issues which lead to fine detail and can result in [[Image noise|noisy]] images even with [[Multisample anti-aliasing|samples]] in the hundreds or thousands.{{Citation needed|date=September 2024}}

Python code:
[[File:Derbail.png|thumb|255x255px|Derbail used on a Julia set of the Burning Ship fractal]]
<!-- Fixed up original version. I then converted to complex in final draft
def pixel(x0: float, y0: float, limit: int=1024) -> int:
    def abs_square(a, b):
        return a * a + b * b

    dbail = 1e6

    x = x0
    y = y0

    dx = 1
    dy = 0

    dx_sum = 0
    dy_sum = 0

    for n in range(1, limit):
        xtemp = x * x - y * y + x0
        y = 2 * x * y + y0
        x = xtemp

        dxtemp = (dx * x - dy * y) * 2 + 1 
        dy = (dy * x + dx * y) * 2 
        dx = dxtemp

        dx_sum += dx
        dy_sum += dy

        if abs_square(dx_sum, dy_sum) >= dbail:
            return n

    return 0
 -->
<syntaxhighlight lang="python3" line="1" start="1">
def mand_der(c0: complex, limit: int=1024):
    def abs_square(c: complex):
        return c.real**2 + c.imag**2

    dbail = 1e6

    c = c0
    dc = 1 + 0j
    dc_sum = 0 + 0j

    for n in range(1, limit):
        c = c**2 + c0
        dc= 2*dc*c + 1
        dc_sum += dc

        if abs_square(dc_sum) >= dbail:
            return n

    return 0
</syntaxhighlight>

==Coloring algorithms==

In addition to plotting the set, a variety of algorithms have been developed to 
* efficiently color the set in an aesthetically pleasing way
* show structures of the data (scientific visualisation)

===Histogram coloring===
<!-- {{Unreferenced section|date=October 2018}} -->
{{More citations needed section|talk="Histogram coloring" section|date=June 2019}}
A more complex coloring method involves using a [[histogram]] which pairs each pixel with said pixel's maximum iteration count before escape/bailout<!--, from 1 to ''n''-->. This method will equally distribute colors to the same overall area, and, importantly, is independent of the maximum number of iterations chosen.<ref>{{cite web |url=http://www.fractalforums.com/programming/newbie-how-to-map-colors-in-the-mandelbrot-set/msg3465/#msg3465 |title=Newbie: How to map colors in the Mandelbrot set? |author= |date=May 2007 |website=www.fractalforums.com |publisher= |access-date=February 11, 2020 |quote= |archive-date=9 September 2022 |archive-url=https://web.archive.org/web/20220909215450/http://www.fractalforums.com/programming/newbie-how-to-map-colors-in-the-mandelbrot-set/msg3465/#msg3465 |url-status=live }}</ref>

This algorithm has four passes. The first pass involves calculating the iteration counts associated with each pixel (but without any pixels being plotted). These are stored in an array: IterationCounts[x][y], where x and y are the x and y coordinates of said pixel on the screen respectively.

{| class="floatright" style="width: 30%; border: 1px solid #c8ccd1; text-align: center; font-size:85%"
|-
|[[File:Mandelbrot-no-histogram-coloring-10000-iterations.png|upright=0.65|center|frameless]]
|[[File:Mandelbrot-no-histogram-coloring-1000-iterations.png|upright=0.65|center|frameless]]
|[[File:Mandelbrot-no-histogram-coloring-100-iterations.png|upright=0.65|center|frameless]]
|-
|[[File:Mandelbrot-histogram-10000-iterations.png|upright=0.65|center|frameless]]
|[[File:Mandelbrot-histogram-1000-iterations.png|upright=0.65|center|frameless]]
|[[File:Mandelbrot-histogram-100-iterations.png|upright=0.65|center|frameless]]
|-
| rowspan=2 colspan=3 style="text-align: center;" |The top row is a series of plots using the escape time algorithm for 10000, 1000 and 100 maximum iterations per pixel respectively. The bottom row uses the same maximum iteration values but utilizes the histogram coloring method. Notice how little the coloring changes per different maximum iteration counts for the histogram coloring method plots.
|}

The first step of the second pass is to create an array of size ''n'', which is the maximum iteration count: NumIterationsPerPixel<!--[''i'']-->. Next, one must iterate over the array of pixel-iteration count pairs, IterationCounts[][], and retrieve each pixel's saved iteration count, ''i'', via e.g. ''i'' = IterationCounts[x][y]. After each pixel's iteration count ''i'' is retrieved, it is necessary to index the NumIterationsPerPixel by ''i'' and increment the indexed value (which is initially zero) -- e.g. NumIterationsPerPixel[''i''] = NumIterationsPerPixel[''i''] + 1 .<!-- This creates the histogram during computation of the image. -->

 '''for''' (x = 0; x < width; x++) '''do'''
     '''for''' (y = 0; y < height; y++) '''do'''
         i:= IterationCounts[x][y]
         NumIterationsPerPixel[i]++

<!--The third pass "renders" each pixel by using the accumulated histogram data--><!-- (Also explained further in the following code sample)-->The third pass iterates through the NumIterationsPerPixel array and adds up all the stored values, saving them in ''total''. The array index represents the number of pixels that reached that iteration count before bailout. <!-- After this step has been completed, the ''total'' value is stored, which is used later. --><!-- Next, each pixel is consulted to determine its maximum iteration count.--><!-- as shown in the following code. -->

 total: = 0
 '''for''' (i = 0; i < max_iterations; i++) '''do'''
     total += NumIterationsPerPixel[i]

After this, the fourth pass begins and all the values in the IterationCounts array are indexed, and, for each iteration count ''i'', associated with each pixel, the count is added to a global sum of all the iteration counts from 1 to ''i'' in the NumIterationsPerPixel array <!--and summed, from 1 to ''i'', up to and including the maximum iteration count of the associated pixel, ''i''--><!--, as shown in the following code-->. This value is then normalized by dividing the sum by the ''total'' value computed earlier.

 hue[][]:= 0.0
 '''for''' (x = 0; x < width; x++) '''do'''
     '''for''' (y = 0; y < height; y++) '''do'''
         iteration:= IterationCounts[x][y]
         '''for''' (i = 0; i <= iteration; i++) '''do'''
             hue[x][y] += NumIterationsPerPixel[i] / total ''/* Must be floating-point division. */''
 
 ...
 
 color = palette[hue[m, n]]
 
 ...

Finally, the computed value is used, e.g. as an index to a color palette.

<!-- If one has a continuous color palette ranging from 0 to 1 for each color component, then the normalized color of each pixel can be computed as follows, using the variables from above. -->
<!-- I don't know who wrote the above, but the "computed as follows" section is missing, so I'm commenting this section out until someone can add the code -->
This method may be combined with the smooth coloring method below for more aesthetically pleasing images.

===Continuous (smooth) coloring===
{{multiple image| direction = vertical
| width = 180
| image1 = Escape Time Algorithm bands.png
| caption1 = This image was rendered with the escape time algorithm. There are very obvious "bands" of color
| image2 = Normalized Iteration Count Algorithm 1.png
| caption2 = This image was rendered with the normalized iteration count algorithm. The bands of color have been replaced by a smooth gradient. Also, the colors take on the same pattern that would be observed if the escape time algorithm were used.
}}
The escape time algorithm is popular for its simplicity. However, it creates bands of color, which, as a type of [[aliasing]], can detract from an image's aesthetic value. This can be improved using an algorithm known as "normalized iteration count",<ref>{{cite journal|last=García|first=Francisco|author2=Ángel Fernández|author3=Javier Barrallo|author4=Luis Martín|title=Coloring Dynamical Systems in the Complex Plane|url=http://math.unipa.it/~grim/Jbarrallo.PDF|accessdate=2008-01-21|archive-date=30 November 2019|archive-url=https://web.archive.org/web/20191130214645/http://math.unipa.it/%7Egrim/Jbarrallo.PDF|url-status=live}}</ref><ref>{{cite web|url=http://linas.org/art-gallery/escape/escape.html|title=Renormalizing the Mandelbrot Escape|author=Linas Vepstas|access-date=11 February 2020|archive-date=14 February 2020|archive-url=https://web.archive.org/web/20200214073730/http://linas.org/art-gallery/escape/escape.html|url-status=live}}</ref> which provides a smooth transition of colors between iterations. The algorithm associates a real number <math>\nu</math> with each value of ''z'' by using the connection of the iteration number with the [[Julia set|potential function]]. This function is given by

: <math>\phi(z) = \lim_{n \to \infty} (\log|z_n|/P^{n}),</math>

where ''z''<sub>''n''</sub> is the value after ''n'' iterations and ''P'' is the power for which ''z'' is raised to in the Mandelbrot set equation (''z''<sub>''n''+1</sub>&nbsp;=&nbsp;''z''<sub>''n''</sub><sup>''P''</sup>&nbsp;+&nbsp;''c'', ''P'' is generally 2).

If we choose a large bailout radius ''N'' (e.g., 10<sup>100</sup>), we have that

: <math>\log|z_n|/P^{n} = \log(N)/P^{\nu(z)}</math>

for some real number <math>\nu(z)</math>, and this is

: <math>\nu(z) = n - \log_P (\log|z_n|/\log(N)),</math>

and as ''n'' is the first iteration number such that |''z''<sub>''n''</sub>| > ''N'', the number we subtract from ''n'' is in the interval <nowiki>[0,&nbsp;1)</nowiki>.

For the coloring we must have a cyclic scale of colors (constructed mathematically, for instance) and containing ''H'' colors numbered from 0 to ''H''&nbsp;−&nbsp;1 (''H'' =&nbsp;500, for instance). We multiply the real number <math>\nu(z)</math> by a fixed real number determining the density of the colors in the picture, take the integral part of this number modulo&nbsp;''H'', and use it to look up the corresponding color in the color table.

For example, modifying the above pseudocode and also using the concept of [[linear interpolation]] would yield
 '''for each''' pixel (Px, Py) on the screen '''do'''
     x0:= scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.5, 1))
     y0:= scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1, 1))
     x:= 0.0
     y:= 0.0
     iteration:= 0
     max_iteration:= 1000
     ''// Here N = 2^8 is chosen as a reasonable bailout radius.''
 
     '''while''' x*x + y*y ≤ (1 << 16) '''and''' iteration < max_iteration '''do'''
         xtemp:= x*x - y*y + x0
         y:= 2*x*y + y0
         x:= xtemp
         iteration:= iteration + 1
 
     ''// Used to avoid floating point issues with points inside the set.''
     '''if''' iteration < max_iteration '''then'''
         ''// sqrt of inner term removed using log simplification rules.''
         log_zn:= log(x*x + y*y) / 2
         nu:= log(log_zn / log(2)) / log(2)
         ''// Rearranging the potential function.''
         ''// Dividing log_zn by log(2) instead of log(N = 1<<8)''
         ''// because we want the entire palette to range from the''
         ''// center to radius 2, NOT our bailout radius.''
         iteration:= iteration + 1 - nu
 
     color1:= palette[floor(iteration)]
     color2:= palette[floor(iteration) + 1]
     ''// iteration % 1 = fractional part of iteration.''
     color:= linear_interpolate(color1, color2, iteration % 1)
     plot(Px, Py, color)

<!--=Further optimizations=-->

=== Exponentially mapped and cyclic iterations ===
[[File:CyclicColoringLch.png|thumb|391x391px|Exponential Cyclic Coloring in LCH color space with shading]]
Typically when we render a fractal, the range of where colors from a given palette appear along the fractal is static. If we desire to offset the location from the border of the fractal, or adjust their palette to cycle in a specific way, there are a few simple changes we can make when taking the final iteration count before passing it along to choose an item from our palette.

When we have obtained the iteration count, we can make the range of colors non-linear. Raising a value normalized to the range [0,1] to a power ''n'', maps a linear range to an exponential range, which in our case can nudge the appearance of colors along the outside of the fractal, and allow us to bring out other colors, or push in the entire palette closer to the border.

<math>v = ((\mathbf{i} / max_i)^\mathbf{S}\mathbf{N})^{1.5} \bmod \mathbf{N}
</math>

where '''i''' is our iteration count after bailout, '''max_i''' is our iteration limit, '''S''' is the exponent we are raising iters to, and '''N''' is the number of items in our palette. This scales the iter count non-linearly and scales the palette to cycle approximately proportionally to the zoom.

We can then plug v into whatever algorithm we desire for generating a color.

=== Passing iterations into a color directly ===
[[File:LCH COLORING.png|thumb|320x320px|Example of exponentially mapped cyclic LCH coloring without shading]]
One thing we may want to consider is avoiding having to deal with a palette or color blending at all. There are actually a handful of methods we can leverage to generate smooth, consistent coloring by constructing the color on the spot.

==== v refers to a normalized exponentially mapped cyclic iter count ====

==== f(v) refers to the [[SRGB#Transfer function ("gamma")|sRGB transfer function]] ====
A naive method for generating a color in this way is by directly scaling '''v''' to 255 and passing it into RGB as such
 rgb = [v * 255, v * 255, v * 255]
One flaw with this is that RGB is non-linear due to gamma; consider linear sRGB instead. Going from RGB to sRGB uses an inverse companding function on the channels. This makes the gamma linear, and allows us to properly sum the colors for sampling.
 srgb = [v * 255, v * 255, v * 255]
[[File:HSV Gradient Example.png|left|thumb|320x320px|HSV Gradient]]

=== HSV coloring ===
HSV Coloring can be accomplished by mapping iter count from [0,max_iter) to [0,360), taking it to the power of 1.5, and then modulo 360.
[[File:HSV Hue Calculation.png|frameless|220x220px]]
We can then simply take the exponentially mapped iter count into the value and return
 hsv = [powf((i / max) * 360, 1.5) % 360, 100, (i / max) * 100]
This method applies to HSL as well, except we pass a saturation of 50% instead.
 hsl = [powf((i / max) * 360, 1.5) % 360, 50, (i / max) * 100]
[[File:LCH Gradient Example.png|left|thumb|320x320px|LCH Gradient]]

=== LCH coloring ===
One of the most perceptually uniform coloring methods involves passing in the processed iter count into LCH. If we utilize the exponentially mapped and cyclic method above, we can take the result of that into the Luma and Chroma channels. We can also exponentially map the iter count and scale it to 360, and pass this modulo 360 into the hue.

<math display="inline">\begin{array}{lcl} 
x & \in & \mathbb{Q+} \\
s_i & =  & (i/max_i)^\mathbf{x} \\
v & = & 1.0 - cos^2(\pi s_i)\\
L & = & 75 - (75v) \\ 
C & = & 28 + (75 - 75v) \\
H & = & (360s_i)^{1.5} \bmod 360
\end{array}
</math>

One issue we wish to avoid here is out-of-gamut colors. This can be achieved with a little trick based on the change in in-gamut colors relative to luma and chroma. As we increase luma, we need to decrease chroma to stay within gamut.
 s = iters/max_i;
 v = 1.0 - powf(cos(pi * s), 2.0);
 LCH = [75 - (75 * v), 28 + (75 - (75 * v)), powf(360 * s, 1.5) % 360];

==Advanced plotting algorithms==
In addition to the simple and slow escape time algorithms already discussed, there are many other more advanced algorithms that can be used to speed up the plotting process.

===Distance estimates===
One can compute the [[distance]] from point ''c'' (in [[exterior (topology)|exterior]] or [[interior (topology)|interior]]) to nearest point on the [[boundary (topology)|boundary]] of the Mandelbrot set.<ref name="Albert Lobo">{{cite web | url=http://www.albertlobo.com/fractals/interior-exterior-distance-bounds-mandelbrot-set | title=Interior and exterior distance bounds for the Mandelbrot set | author=Albert Lobo | access-date=29 April 2021 | archive-date=9 September 2022 | archive-url=https://web.archive.org/web/20220909215805/http://www.albertlobo.com/fractals/interior-exterior-distance-bounds-mandelbrot-set | url-status=live }}</ref><ref>{{cite web |url=http://www.imajeenyus.com/mathematics/20121112_distance_estimates/distance_estimation_method_for_fractals.pdf |title=Distance estimation method for drawing Mandelbrot and Julia sets |last=Wilson |first=Dr. Lindsay Robert |date=2012 |access-date=3 May 2021 |archive-date=3 May 2021 |archive-url=https://web.archive.org/web/20210503103652/http://www.imajeenyus.com/mathematics/20121112_distance_estimates/distance_estimation_method_for_fractals.pdf |url-status=live }}</ref>

====Exterior distance estimation====
The proof of the [[Mandelbrot set#Basic properties|connectedness]] of the Mandelbrot set in fact gives a formula for the [[External ray#Uniformization|uniformizing map]] of the [[complement (set theory)|complement]] of <math>M</math> (and the [[derivative]] of this map). By the [[Koebe quarter theorem]], one can then estimate the distance between the midpoint of our [[pixel]] and the Mandelbrot set up to a factor of 4.

In other words, provided that the maximal number of iterations is sufficiently high, one obtains a picture of the Mandelbrot set with the following properties:

# Every pixel that contains a point of the Mandelbrot set is colored black.
# Every pixel that is colored black is close to the Mandelbrot set.

[[File:Demm 2000 Mandelbrot set.jpg|thumb|Exterior distance estimate may be used to color whole complement of Mandelbrot set]]
The upper bound ''b'' for the distance estimate of a pixel ''c'' (a complex number) from the Mandelbrot set is given by<ref name="Cheritat">{{cite web |url=https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Boundary_detection_methods_via_distance_estimators |title=Boundary detection methods via distance estimators |last=Chéritat |first=Arnaud |author-link=Arnaud Chéritat |date=2016 |access-date=2 January 2023 |archive-date=18 December 2022 |archive-url=https://web.archive.org/web/20221218100822/https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set |url-status=live }}</ref><ref name="Syntopia">{{cite web |url=http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations |title=Distance Estimated 3D Fractals (V): The Mandelbulb & Different DE Approximations |last=Christensen |first=Mikael Hvidtfeldt |date=2011 |access-date=10 May 2021 |archive-date=13 May 2021 |archive-url=https://web.archive.org/web/20210513033347/http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations/ |url-status=live }}</ref><ref>{{cite book |title=Hypercomplex Iterations: Distance Estimation and Higher Dimensional Fractals |last=Dang |first=Yumei |author2=Louis Kauffman |author-link2=Louis Kauffman |author3=Daniel Sandin |author-link3=Daniel J. Sandin |date=2002 |pages=17–18 |chapter=Chapter 3.3: The Distance Estimation Formula |publisher=World Scientific |url=https://www.evl.uic.edu/hypercomplex/html/book/book.pdf |access-date=29 April 2021 |archive-date=23 March 2021 |archive-url=https://web.archive.org/web/20210323153331/https://www.evl.uic.edu/hypercomplex/html/book/book.pdf |url-status=live }}</ref>

: <math>b=\lim_{n \to \infty} \frac{2 \cdot |{P_c^n(c)| \cdot \ln|{P_c^n(c)}}|}{|\frac{\partial}{\partial{c}} P_c^n(c)|},</math>
where 
* <math>P_c(z) \,</math> stands for [[complex quadratic polynomial]]
* <math>P_c^n(c)</math> stands for ''n'' iterations of <math>P_c(z) \to z</math> or <math>z^2 + c \to z</math>, starting with <math>z=c</math>:<!--
 --> <math>P_c^{ 0}(c) = c</math>,<!--
 --> <math>P_c^{ n+1}(c) = P_c^n(c)^2 + c</math>;
* <math>\frac{\partial}{\partial{c}} P_c^n(c)</math> is the derivative of<!--
 --> <math>P_c^n(c)</math> with respect to c. This derivative can be found by starting with<!--
 --> <math>\frac{\partial}{\partial{c}} P_c^{ 0}(c) = 1</math><!--
 --> and then <math>\frac{\partial}{\partial{c}} P_c^{ n+1}(c) = 2\cdot{}P_c^n(c)\cdot\frac{\partial}{\partial{c}} P_c^n(c) + 1</math>.<!--
 --> This can easily be verified by using the chain rule for the derivative.

The idea behind this formula is simple: When the [[equipotential]] lines for the potential function <math>\phi(z)</math> lie close, the number <math>|\phi'(z)|</math> is large, and conversely, therefore the equipotential lines for the function <math>\phi(z)/|\phi'(z)|</math> should lie approximately regularly.

From a mathematician's point of view, this formula only works in limit where ''n'' goes to infinity, but very reasonable estimates can be found with just a few additional iterations after the main loop exits.

Once ''b'' is found, by the Koebe 1/4-theorem, we know that there is no point of the Mandelbrot set with distance from ''c'' smaller than ''b/4''.

The distance estimation can be used for drawing of the boundary of the Mandelbrot set, see the article [[Julia set]]. In this approach, pixels that are sufficiently close to M are drawn using a different color. This creates drawings where the thin "filaments" of the Mandelbrot set can be easily seen. This technique is used to good effect in the B&W images of Mandelbrot sets in the books "The Beauty of Fractals<ref>{{cite book |title=The Beauty of Fractals |last=Peitgen |first=Heinz-Otto |author2=Richter Peter |year=1986 |publisher=Springer-Verlag |location=Heidelberg |isbn=0-387-15851-0 |pages= |title-link=The Beauty of Fractals }}</ref>" and "The Science of Fractal Images".<ref>{{cite book |title=The Science of Fractal Images |last=Peitgen |first=Heinz-Otto |author2=Saupe Dietmar |year=1988 |publisher=Springer-Verlag |location=New York |isbn=0-387-96608-0 |pages=202 |title-link=The Beauty of Fractals }}</ref>

Here is a sample B&W image rendered using Distance Estimates:
[[File:Mandel zoom 06 double hook B&W DE.jpg|thumb|none|This is a B&W image of a portion of the Mandelbrot set rendered using Distance Estimates (DE)]]

Distance Estimation can also be used to render [[Mandelbrot Set#3D images of Mandelbrot and Julia sets|3D images of Mandelbrot and Julia sets]]

====Interior distance estimation====
[[File:Mandelbrot Interior 600.png|thumb|Pixels colored according to the estimated interior distance]]
It is also possible to estimate the distance of a limitly periodic (i.e., [[Mandelbrot set#Hyperbolic components|hyperbolic]]) point to the boundary of the Mandelbrot set. The upper bound ''b'' for the distance estimate is given by<ref name="Albert Lobo" />

: <math>b=\frac{1-\left|{\frac{\partial}{\partial{z}}P_c^p(z_0)}\right|^2}
  {\left|{\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}P_c^p(z_0) +
        \frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}P_c^p(z_0)
        \frac{\frac{\partial}{\partial{c}}P_c^p(z_0)}
             {1-\frac{\partial}{\partial{z}}P_c^p(z_0)}} \right|},</math>
where
* <math>p</math> is the period,
* <math>c</math> is the point to be estimated,
* <math>P_c(z)</math> is the [[complex quadratic polynomial]] <math>P_c(z)=z^2 + c</math>
* <math>P_c^p(z_0)</math> is the <math>p</math>-fold iteration of <math>P_c(z) \to z</math>, starting with <math>P_c^{ 0}(z) = z_0</math>
* <math>z_0</math> is any of the <math>p</math> points that make the [[attractor]] of the iterations of <math>P_c(z) \to z</math> starting with<!--
 --> <math>P_c^{ 0}(z) = c</math>; <math>z_0</math> satisfies <math>z_0 = P_c^p(z_0)</math>,
* <math>\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}P_c^p(z_0)</math>,<!--
 --> <math>\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}P_c^p(z_0)</math>,<!--
 --> <math>\frac{\partial}{\partial{c}}P_c^p(z_0)</math> and<!--
 --> <math>\frac{\partial}{\partial{z}}P_c^p(z_0)</math> are various derivatives of<!--
 --> <math>P_c^p(z)</math>, evaluated at <math>z_0</math>.

Analogous to the exterior case, once ''b'' is found, we know that all points within the distance of ''b''/4 from ''c'' are inside the Mandelbrot set.

There are two practical problems with the interior distance estimate: first, we need to find <math>z_0</math> precisely, and second, we need to find <math>p</math> precisely.
The problem with <math>z_0</math> is that the convergence to <math>z_0</math> by iterating <math>P_c(z)</math> requires, theoretically, an infinite number of operations.
The problem with any given <math>p</math> is that, sometimes, due to rounding errors, a period is falsely identified to be an integer multiple of the real period (e.g., a period of 86 is detected, while the real period is only 43=86/2). In such case, the distance is overestimated, i.e., the reported radius could contain points outside the Mandelbrot set.
[[File:MandelbrotOrbitInfimum.png|thumb|3D view: smallest absolute value of the orbit of the interior points of the Mandelbrot set]]

===Cardioid / bulb checking===
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 that:

:<math> p = \sqrt{  \left(x - \frac{1}{4}\right)^2 + y^2} </math>,
:<math> x \leq p - 2p^2 + \frac{1}{4} </math>,
:<math> (x+1)^2 + y^2 \leq \frac{1}{16} </math>,

where ''x'' represents the real value of the point and ''y'' the imaginary value. The first two equations determine that the point is within the cardioid, the last the period-2 bulb.

The cardioid test can equivalently be performed without the square root:

:<math> q = \left(x - \frac{1}{4}\right)^2 + y^2, </math>
:<math> q \left(q + \left(x - \frac{1}{4}\right)\right) \leq \frac{1}{4}y^2. </math>

3rd- and higher-order buds do not have equivalent tests, because they are not perfectly circular.<ref>{{cite web
|url=http://linas.org/art-gallery/bud/bud.html
|title=Mandelbrot Bud Maths
|access-date=11 February 2020
|archive-date=14 February 2020
|archive-url=https://web.archive.org/web/20200214073800/http://linas.org/art-gallery/bud/bud.html
|url-status=live
}}</ref> However, it is possible to find whether the points are within circles inscribed within these higher-order bulbs, preventing many, though not all, of the points in the bulb from being iterated.

=== Periodicity checking ===

To prevent having to do huge numbers of iterations for points inside the set, one can perform periodicity checking, which checks whether a point reached in iterating a pixel has been reached before. If so, the pixel cannot diverge and must be in the set. Periodicity checking is a trade-off, as the need to remember points costs data management instructions and memory, but saves computational instructions. However, checking against only one previous iteration can detect many periods with little performance overhead. For example, within the while loop of the pseudocode above, make the following modifications:

 xold := 0
 yold := 0
 period := 0
 '''while''' (x*x + y*y ≤ 2*2 '''and''' iteration < max_iteration) '''do'''
     xtemp := x*x - y*y + x0
     y := 2*x*y + y0
     x := xtemp
     iteration := iteration + 1 
  
     '''if''' x ≈ xold '''and''' y ≈ yold '''then'''
         iteration := max_iteration  /* Set to max for the color plotting */
         break        /* We are inside the Mandelbrot set, leave the while loop */
  
     period:= period + 1
     '''if''' period > 20 '''then'''
         period := 0
         xold := x
         yold := y

The above code stores away a new x and y value on every 20th iteration, thus it can detect periods that are up to 20 points long.

===Border tracing / edge checking===
[[File:Mandelbrot DEM Sobel.png|thumbnail|right|Edge detection using [[Sobel filter]] of hyperbolic components of Mandelbrot set]]Because the Mandelbrot set is [[Full set (topology)|full]],<ref name="Proof of connectedness">{{cite journal |last1=Douady |first1=Adrien |last2=Hubbard |first2=John |title=Exploring the Mandelbrot set. Exploring the Mandelbrot set. The Orsay Notes. |date=2009 |url=https://www.researchgate.net/publication/228540477 |access-date=9 April 2023}}</ref> any point enclosed by a closed shape whose borders lie entirely within the Mandelbrot set must itself be in the Mandelbrot set. Border tracing works by following the [[polynomial lemniscate|lemniscates]] of the various iteration levels (colored bands) all around the set, and then filling the entire band at once. This also provides a speed increase because large numbers of points can be now skipped.<ref>{{cite web|url=http://www.reocities.com/CapeCanaveral/5003/mandel.htm|title=Boundary Tracing Method|archiveurl=https://web.archive.org/web/20150220012221/http://www.reocities.com/CapeCanaveral/5003/mandel.htm|archivedate=2015-02-20}}</ref>

[[File:Boundary following movie.webm|thumb|Animation of border tracing]]

In the animation shown, points outside the set are colored with a 1000-iteration escape time algorithm. Tracing the set border and filling it, rather than iterating the interior points, reduces the total number of iterations by 93.16%. With a higher iteration limit the benefit would be even greater.

=== Rectangle checking ===

Rectangle checking is an older and simpler method for plotting the Mandelbrot set. The basic idea of rectangle checking is that if every pixel in a rectangle's border shares the same amount of iterations, then the rectangle can be safely filled using that number of iterations. There are several variations of the rectangle checking method, however, all of them are slower than the border tracing method because they end up calculating more pixels. One variant just calculates the corner pixels of each rectangle, however, this causes damaged pictures more often than calculating the entire border, thus it only works reasonably well if only small boxes of around 6x6 pixels are used, and no recursing in from bigger boxes. ([[Fractint]] method.)

The most simple rectangle checking method lies in checking the borders of equally sized rectangles, resembling a grid pattern. (Mariani's algorithm.)<ref>{{cite magazine|title=Computer Recreations, February 1989; A tour of the Mandelbrot set aboard the Mandelbus |last=Dewdney |first=A. K. |year=1989 |magazine=Scientific American |pages=111|jstor=24987149 }} {{subscription required}}</ref>

A faster and slightly more advanced variant is to first calculate a bigger box, say 25x25 pixels. If the entire box border has the same color, then just fill the box with the same color. If not, then split the box into four boxes of 13x13 pixels, reusing the already calculated pixels as outer border, and sharing the inner "cross" pixels between the inner boxes. Again, fill in those boxes that has only one border color. And split those boxes that don't, now into four 7x7 pixel boxes. And then those that "fail" into 4x4 boxes. (Mariani-Silver algorithm.)

Even faster is to split the boxes in half instead of into four boxes. Then it might be optimal to use boxes with a 1.4:1 [[aspect ratio]], so they can be split like [[ISO 216|how A3 papers are folded]] into A4 and A5 papers. (The DIN approach.)

As with border tracing, rectangle checking only works on areas with one discrete color. But even if the outer area uses smooth/continuous coloring then  rectangle checking will still speed up the costly inner area of the Mandelbrot set. Unless the inner area also uses some smooth coloring method, for instance [[#Interior distance estimation|interior distance estimation]].

=== Symmetry utilization ===
The horizontal symmetry of the Mandelbrot set allows for portions of the rendering process to be skipped upon the presence of the real axis in the final image. However, regardless of the portion that gets mirrored, the same number of points will be rendered.

Julia sets have symmetry around the origin. This means that quadrant 1 and quadrant 3 are symmetric, and quadrants 2 and quadrant 4 are symmetric. Supporting symmetry for both Mandelbrot and Julia sets requires handling symmetry differently for the two different types of graphs.

=== Multithreading ===

Escape-time rendering of Mandelbrot and Julia sets lends itself extremely well to parallel processing. On multi-core machines the area to be plotted can be divided into a series of rectangular areas which can then be provided as a set of tasks to be rendered by a pool of rendering threads. This is an [[embarrassingly parallel]]<ref>http://courses.cecs.anu.edu.au/courses/COMP4300/lectures/embParallel.4u.pdf {{Webarchive|url=https://web.archive.org/web/20200127230256/http://courses.cecs.anu.edu.au/courses/COMP4300/lectures/embParallel.4u.pdf |date=27 January 2020 }} {{Bare URL PDF|date=March 2022}}</ref> computing problem. (Note that one gets the best speed-up by first excluding symmetric areas of the plot, and then dividing the remaining unique regions into rectangular areas.)<ref>http://cseweb.ucsd.edu/groups/csag/html/teaching/cse160s05/lectures/Lecture14.pdf {{Webarchive|url=https://web.archive.org/web/20200126040920/http://cseweb.ucsd.edu/groups/csag/html/teaching/cse160s05/lectures/Lecture14.pdf |date=26 January 2020 }} {{Bare URL PDF|date=March 2022}}</ref>

Here is a short video showing the Mandelbrot set being rendered using multithreading and symmetry, but without boundary following:

[[File:Multi-threaded mandelbrot render w symmetry.webm|thumb|none|This is a short video showing rendering of a Mandelbrot set using multi-threading and symmetry, but with boundary following turned off.]]

Finally, here is a video showing the same Mandelbrot set image being rendered using multithreading, symmetry, '''and''' boundary following:

[[File:Multi-threaded mandelbrot render w bounary following & symmetry.webm|thumb|none|This is a short video showing rendering of a Mandelbrot set using boundary following, multi-threading, and symmetry]]

<!--== Advanced bailout method ==
Simple programs and scripts generally tend to set the escape value to two. This process can be improved by utilizing the distance from the origin and the point being rendered with the Pythagorean theorem by summing the squares of the real and imaginary portions of <math>z</math> and then escaping if the value is larger or equal to four. The result of this optimization is a faster rendering of the image.-->

===Perturbation theory and series approximation===
Very highly magnified images require more than the standard 64–128 or so bits of precision that most hardware [[floating-point unit]]s provide, requiring renderers to use slow "BigNum" or "[[Arbitrary-precision arithmetic|arbitrary-precision]]" math libraries to calculate. However, this can be sped up by the exploitation of [[perturbation theory]]. Given

:<math> z_{n+1} = z_n^2 + c </math>

as the iteration, and a small epsilon and delta, it is the case that

:<math> (z_n + \epsilon)^2 + (c + \delta) = z_n^2 + 2z_n\epsilon + \epsilon^2 + c + \delta, </math>

or

:<math> = z_{n+1} + 2z_n\epsilon + \epsilon^2 + \delta, </math>

so if one defines

:<math> \epsilon_{n+1} = 2z_n\epsilon_n + \epsilon_n^2 + \delta, </math>

one can calculate a single point (e.g. the center of an image) using high-precision arithmetic (''z''), giving a ''reference orbit'', and then compute many points around it in terms of various initial offsets delta plus the above iteration for epsilon, where epsilon-zero is set to 0. For most iterations, epsilon does not need more than 16 significant figures, and consequently hardware floating-point may be used to get a mostly accurate image.<ref>{{cite web
|url=http://www.fractalforums.com/announcements-and-news/superfractalthing-arbitrary-precision-mandelbrot-set-rendering-in-java/
|title=Superfractalthing - Arbitrary Precision Mandelbrot Set Rendering in Java
|access-date=11 February 2020
|archive-date=30 June 2020
|archive-url=https://web.archive.org/web/20200630043303/http://www.fractalforums.com/announcements-and-news/superfractalthing-arbitrary-precision-mandelbrot-set-rendering-in-java/
|url-status=live
}}</ref> There will often be some areas where the orbits of points diverge enough from the reference orbit that extra precision is needed on those points, or else additional local high-precision-calculated reference orbits are needed. By measuring the orbit distance between the reference point and the point calculated with low precision, it can be detected that it is not possible to calculate the point correctly, and the calculation can be stopped. These incorrect points can later be re-calculated e.g. from another closer reference point.

Further, it is possible to approximate the starting values for the low-precision points with a truncated [[Taylor series]], which often enables a significant amount of iterations to be skipped.<ref>{{cite journal
|author=K. I. Martin
|url=http://www.superfractalthing.co.nf/sft_maths.pdf
|title=Superfractalthing Maths
|journal=
|access-date=11 February 2020
|archive-url=https://web.archive.org/web/20140628114658/http://www.superfractalthing.co.nf/sft_maths.pdf
|archive-date=28 June 2014
|url-status=dead
}}</ref>
Renderers implementing these techniques are [[Kalles Fraktaler|publicly available]] and offer speedups for highly magnified images by around two orders of magnitude.<ref>{{cite web
|url=http://www.chillheimer.de/kallesfraktaler/
|title=Kalles Fraktaler 2
|access-date=11 February 2020
|archive-date=24 February 2020
|archive-url=https://web.archive.org/web/20200224054206/http://www.chillheimer.de/kallesfraktaler/
|url-status=live
}}</ref>

An alternate explanation of the above:

For the central point in the disc <math> c </math> and its iterations <math> z_n </math>, and an arbitrary point in the disc <math> c + \delta </math> and its iterations <math> z'_n </math>, it is possible to define the following iterative relationship:

:<math> z'_{n} = z_{n} + \epsilon_{n} </math>

With <math> \epsilon_{1} = \delta </math>. Successive iterations of <math> \epsilon_n </math> can be found using the following:

:<math> z'_{n+1} = {z'_n}^2 + (c + \delta) </math>

:<math> z'_{n+1} = (z_n + \epsilon_n)^2 + c + \delta </math>

:<math> z'_{n+1} = {z_n}^2 + c + 2z_n\epsilon_n + {\epsilon_n}^2 + \delta</math>

:<math> z'_{n+1} = z_{n+1} + 2z_n\epsilon_n + {\epsilon_n}^2 + \delta</math>

Now from the original definition:
:<math> z'_{n+1} = z_{n+1} + \epsilon_{n+1} </math>,

It follows that:

:<math> \epsilon_{n+1} = 2z_n\epsilon_n + {\epsilon_n}^2 + \delta </math>

As the iterative relationship relates an arbitrary point to the central point by a very small change <math> \delta </math>, then most of the iterations of <math> \epsilon_n </math> are also small and can be calculated using floating point hardware.

However, for every arbitrary point in the disc it is possible to calculate a value for a given <math> \epsilon_{n} </math> without having to iterate through the sequence from <math> \epsilon_0 </math>, by expressing <math> \epsilon_n </math> as a [[power series]] of <math> \delta </math>.

:<math> \epsilon_n = A_{n}\delta + B_{n}\delta^2 + C_{n}\delta^3 + \dotsc </math>

With <math> A_{1} = 1, B_{1} = 0, C_{1} = 0, \dotsc </math>.

Now given the iteration equation of <math> \epsilon </math>, it is possible to calculate the coefficients of the power series for each <math> \epsilon_n </math>:

:<math> \epsilon_{n+1} = 2z_n\epsilon_n + {\epsilon_n}^2 + \delta </math>

:<math> \epsilon_{n+1} = 2z_n(A_n\delta + B_n\delta^2 + C_n\delta^3 + \dotsc) + (A_n\delta + B_n\delta^2 + C_n\delta^3 + \dotsc)^2 + \delta </math>

:<math> \epsilon_{n+1} = (2z_nA_n+1)\delta + (2z_nB_n + {A_n}^2)\delta^2 + (2z_nC_n + 2A_nB_n)\delta^3  + \dotsc </math>

Therefore, it follows that:

:<math> A_{n+1} = 2z_nA_n + 1 </math>
:<math> B_{n+1} = 2z_nB_n + {A_n}^2 </math>
:<math> C_{n+1} = 2z_nC_n + 2A_nB_n </math>
:<math> \vdots </math>

The coefficients in the power series can be calculated as iterative series using only values from the central point's iterations <math> z </math>, and do not change for any arbitrary point in the disc. If <math> \delta </math> is very small, <math> \epsilon_n </math> should be calculable to sufficient accuracy using only a few terms of the power series. As the Mandelbrot Escape Contours are 'continuous' over the complex plane, if a points escape time has been calculated, then the escape time of that points neighbours should be similar. Interpolation of the neighbouring points should provide a good estimation of where to start in the <math> \epsilon_n </math> series.

Further, separate interpolation of both real axis points and imaginary axis points should provide both an upper and lower bound for the point being calculated. If both results are the same (i.e. both escape or do not escape) then the difference <math> \Delta n </math> can be used to recuse until both an upper and lower bound can be established. If floating point hardware can be used to iterate the <math> \epsilon </math> series, then there exists a relation between how many iterations can be achieved in the time it takes to use BigNum software to compute a given <math> \epsilon_n </math>. If the difference between the bounds is greater than the number of iterations, it is possible to perform binary search using BigNum software, successively halving the gap until it becomes more time efficient to find the escape value using floating point hardware.

==References==
{{Reflist|30em}}

{{Fractals}}
{{Mathematics and art}}
{{Visualization}}
{{Computer science}}
{{Computer graphics}}

{{DEFAULTSORT:Mandelbrot Set}}
[[Category:Fractals]]
[[Category:Articles with example pseudocode]]
[[Category:Complex dynamics]]
[[Category:Graphics software]]
[[Category:Computer graphics]]
[[Category:Algorithms]]
[[Category:Articles containing video clips]]