This tutorial introduces how to make patterns with complex numbers in Blender. It builds off of tutorials on creative coding in Blender and scripting Cycles nodes. For that reason, it assumes some prior acquaintance with Blender’s Python API and Open Shading Language (OSL).
In this tutorial, we’ll go over
- the maths behind complex numbers;
- how to define complex number operations with nodes, Python and OSL;
- how to create 2D patterns (the Mobius transformation, Mandelbrot set, Julia set, Pickover stalks and tricorn).
- how to create a 3D Mandelbulb.
Last, we’ll look briefly at how a fractal pattern can be mixed with other patterns (for shaders) and modifiers (for meshes generated by Python).
Since the details of implementing complex number functions can take a while, readers are encouraged to skip the initial section of this tutorial where they are defined, look over what patterns can be made with these functions, then return if interested.
This tutorial was written for Blender 2.8 beta.
First, a refresher on the math behind complex numbers. A complex number is written as a+bi. a and b are real scalars. i is an imaginary number, the square root of -1. When i is raised to an even whole number -1 results,
with zero being the exception. When graphed with an Argand Diagram, i.e., on the complex plane, the real axis is horizontal and the imaginary axis is vertical. That means a complex number z
0.3248+0.1875i would look like so:
A complex number can also be expressed in polar coordinates. The number z is expressed as a phase of
30 degrees and an absolute of
The angle is called the argument in some implementations, the phase in others. The radius is alternatively called the modulus or absolute (
abs). In math notation, the absolute is signaled with one or two vertical bars.
The phase is often symbolized by the Greek letter phi, φ.
Notes on Implementation
The similarity between complex numbers and two-dimensional (2D) vectors means that vectors can be used to store and to visualize them. In Cartesian coordinates, the x component carries the real number; the y component, the imaginary. A complex number’s argument matches the vector’s heading; the modulus, a vector’s length or magnitude. For vector implementations where 2D vectors are not distinguished from 3D, care must taken that the z component be treated as zero.
We’ll be implementing complex number operations in Python, OSL and Blender Cycles nodes until we get to pattern-making. A few notes on how these languages will influence our implementation:
- As of Blender 2.8, Cycles math nodes include
sqrt. Nodes are not conducive to iterative algorithms, so we won’t be creating fractal patterns with nodes.
- OSL’s language specification remains the best source when in doubt.
structss can be defined in OSL with Blender 2.8. Make sure to convert from and to recognized data types, rather than use a
structas an input or output. For backwards compatibility, we’ll not be using this in our examples, but the approach looks like so:
- Python handles complex numbers out of the box. Just as real number functions are found in the
mathlibrary, so too with complex functions in the
numpyalso supports complex operations, but we will not be using it for this tutorial.
- Blender’s Python API has been updated for 2.8. Operations between data types in mathutils now use an
@instead of an
*operator. For vectors
a @ bis the dot product. For matrix
m @ vapplies the affine transformation
v. For quaternion
q @ vrotates
Those differences noted, we break complex operations down into four stages:
- Unary operations, pt. 1: absolute, phase, reciprocal;
- Conversion between rectilinear and polar coordinates;
- Binary operations, pt. 1: add, subtract, multiply, divide, power;
- Unary operations, pt. 2: exponent, logarithm, cosine, sine.
We adopt this order because conversions between polar and rectilinear (or Cartesian) coordinates depend on the absolute and phase; division depends on the reciprocal; exponent and logarithm depend on polar-rectilinear conversions. Last, the ability to create a fractal pattern depends on exponent and logarithm at work in the complex power operation.
Unary Operations, Pt. 1
In this section, we cover
Abs Sq / Modulus Sq / R Sq
separate xyzNeeded for:
This is analogous to a vector’s magnitude squared. If we know that the input vector
a has a
z component equal to
0.0, then we can use
dot(a, a) in OSL. With
Vectors in Python, we have the expressions
a @ a,
Vector.dot(a, a). This can be used to avoid the
sqrt implicit in
Abs / Modulus / R
mandelbrot (and variants),
This is analogous to a vector’s magnitude. If we know that the input
a has a
z component equal to
0.0, then we can use
length(a) (OSL) or
a.magnitude (Python). In Python,
len(a) returns the number of components a vector has, not its magnitude; for example, a 2D vector has 2 components,
y. With Cycles nodes, if
sqrt is not available,
pow(abssq(z), 0.5) can be used.
pow(dot(a, a), 0.5) will work if we know that
In the context of visualizing fractals, the
abs function is key to finding when the complex number escapes, and thus when the iteration should cease. The usual test is that
z has escaped when
abs(z) > 2.0.
Phase / Arg / Phi
separate xyzNeeded for:
atan2 takes the
imag component as its first argument; the
real component as its second. It returns a value in the range [-π, π]. Use
mod (OSL) or
% (Python) to shift it to [0, τ]. The
Modulo function in the math node is based off
floor, and works more like OSL’s
fmod. For older versions of Blender, we need to use
atan and handle the domain of our inputs as expressed in the formula
(modified from the Wikipedia entry on complex numbers). Shader languages typically handle division-by-zero by returning zero regardless of the numerator’s value. Compare this behavior to throwing an exception (Python) or returning
abs is calculated along the way to finding the
phase, we can include it as an output, making this the same as
polar below. The function depends on two others. The first,
neq (not equals) is defined as
In scenarios where real numbers store Boolean (true or false) values,
maximum can be used for logical
OR. The second,
mix (unclamped linear interpolation) is defined as
mix acts like a ternary operator. For example, the expression
float a = c > 0.0 ? 2.0 : 3.0; is analogous to
float a = mix(3.0, 2.0, c > 0.0);.
mix isn't needed when one of the alternatives is
0.0. This is the case with
float a = mix(0.0, M_PI, real < 0.0); can be simplified to
float a = M_PI * (real < 0.0);.
conjugate is usually too trivial to implement. It can, however, be an extra layer of abstraction to clarify
In math notation, the conjugate of z has a horizontal bar over it. The
conjugate can also be used to find the
as will become clearer after we address complex multiplication.
Reciprocal / Inverse
The reciprocal is most useful in defining division of one complex number by another. Sometimes referred to as an inverse, its math notation is a superscript negative one.
Given that complex division is already defined in Python, we won’t need to use it directly in that environment.
In this section we cover the conversions
Binary Operations, Pt. 1
In this section, we cover
div. We also define
pow, though it depends on
log, which have not yet been defined.
mandelbrot (and variants),
separate xyzNeeded for:
mandelbrot (and variants)
Unary Operations, Pt. 2
This section covers
separate xyzNeeded for:
To implement the above in Cycles nodes, we need hyperbolic cosine (
and hyperbolic sine (
To offer a sense of the pattern created by
cos, we use a diagnostic texture on a plane to animate the following:
Once we start with fractals, we’ll stop using Cycles nodes. This is because fractal algorithms require us to iterate a core algorithm.
Since we know that the vectors returned by
cMult will have
z set to
0.0, we can safely use vector addition.
The animated gif, above, illustrates the type of pattern generated by this algorithm when we plug in an object’s coordinates to z. The manner in which the constants a, b, c and d are changed yields different transformations.
Then, we need a general purpose function that adds an object to a collection (the scene’s main collection by default). This object will be an
Empty if no data is supplied.
The next function we need sets the new object’s transform. Blender’s quaternion takes
w, the real component, as its first argument, not its last. The identity quaternion
(1.0, 0.0, 0.0, 0.0) is equal to Euler angles
(0.0, 0.0, 0.0) and an axis-angle of
0.0about the axis
(0.0, 1.0, 0.0). Quaternions are used because they allow animation without gimbal lock.
Next, we create a lattice data by script. The lattice contains a list of points, each of which is a LatticePoint. The coordinates u, v, and w correspond to x, y and z. We use them when determining the dimensions, or resolution.
We next create a function which loops through the original points in the lattice — which are read-only — and then setting the lattice’s deformed point to the results of the Mobius transformation.
Since the 2D nature of the function nullifies any depth, we can restore some volume by setting the z coordinate by an arbitrary value derived from the x and y dimension.
With a lattice in place, any object with a lattice modifier that targets that lattice will be deformed accordingly. The lattice modifier includes a
Strength factor, ranging from [0, 1] which can be animated with keyframes.
One of the most iconic fractals generated by complex numbers, the Mandelbrot set is named after Benoit Mandelbrot. Mandelbrot’s discovery was guided by a belief that roughness, not only the smooth curves of calculus, is fundamental to natural patterns. In mathematical notation, the function we iterate looks like so:
where c is a complex constant. This function is repeated until |z| exceeds the upper limit,
2.0, or the maximum number of iterations is reached. This can be generalized by replacing the exponent 2 with a variable, labeled e below.
If we want, we can subdivide this variation into smaller ones: where e is a positive integer (
4…), e is a positive real number (
2.2…), e is a signed real number (
1.5…). Paul Bourke’s article “Mandelbrot at higher powers” offers a nice gallery of these variations. We’ll recreate higher-power Mandelbrot sets from there.
floats in place of Booleans; when
ExcludeUpper is equal to zero, it is
ExcludeUpper is true and
i has reached the maximum before escaping, then
Fac is left at
0.0. In this case, the bulb’s interior remains black. The global variable
P is the point currently being sampled.
As Bourke observes, raising the fractal to a positive whole number
b results in
b - 1 bulbs rotated around a pivot. Raising the fractal to a power with a fractional component results in horizontal mirroring. The maximum number of iterations is kin to the ‘refinement’ of the pattern. Set at a low number, the visual is a blob. The OSL script can be hooked up to other nodes, like so:
Fac is calculated by dividing the last iteration by the maximum allowed. This can lead to color banding; we won’t implement a smoothing algorithm, but it is discussed at Wikipedia and by Inigo Quilez. Not all programming languages support
log by an arbitrary base; OSL does. This means code such as
log(a) / log(b); can be condensed to
log(a, b);. A different algorithm, based on the Lyapunov exponent, is needed to depict negative exponents properly; see this discussion for more.
Z at its last iteration and its
Abs are output by this shader as well. The pattern formed by
Z is illustrated below (it has been converted to a
color with the formula
0.5 + 0.5 * normalize(z);).
To visualize this in Python, we’ll create a 2D grid that we populate with cubes.
First we break out the Mandelbrot algorithm into its own function. This function isolates the conversion between complex numbers and vectors. In the function definitions above, we used complex literals; to create a complex number in Python from preexisting numbers, we use
To match an OSL shader’s ability to return multiple outputs, any value with aesthetic potential is included in a dictionary. This function is then called by one which manages the
points is expected to contain vectors. We will create this list from nested for-loops:
The outer loop runs through each row, and therefore determines the y coordinate. The inner loop runs through each column, and determines the x coordinate. The
scale variables refer to the complex plane.
To visualize the fractal with cubes, we work with
data instead of
ops or BMesh functions. For a more detailed walk-through of this approach, Diego Gangl’s “Meshes with Python & Blender” is a handy reference.
First, we create a list of six tuples. The four indices per tuple indicate the corners of a quadrilateral face on the cube. These indices will be used to find a vector in a list that we will be creating later.
Since a mesh data containing multiple cubes will store all the faces in one cumulative list, we allow an offset to be added each time we generate a new set of six face indices. We then generate the eight vertices of a unit cube:
For a cube that scales and rotates from an arbitrary pivot, the function accepts a
pivot which is subtracted from each vector. We use a matrix to transform the vertices. The matrix workflow consists of creating separate matrices from
static functions (Translation, Rotation, Scale), then compositing them via multiplication. We won’t be rotating cubes in this example.
Matrix.Scale doesn’t do nonuniform scaling, so we roll our own:
Matrix is constructed using vectors with 4 components each. With these three functions in place, we create the following to visualize the results of the earlier
More on from_pydata and validate can be found in the Python reference. We depend on
append_obj_from_data, defined earlier in the Mobius transformation section. The easing function which determines the relative size of a cube based on the return factor, by default, is
An alternative example,
float_smooth_step is also included. One last function puts all the above into action:
We have to be careful when scripting due to the high number of iterations required. A Python script could easily tank Blender. For example, a 32 x 32 grid has 1024 loci, which are fed to the
mandelbrot function as
center. Supposing a worst case scenario — where the for-loop reaches its max iteration — a fractal refined to 150 iterations would require 153,600 calls to
mandelbrot_step. One step we can take in
visualize_fractal is to prune away cubes where the factor is out of bounds.
The code above is written for parity with OSL and to illustrate the relationship between complex numbers, fractals and vectors. It’s not written for performance, but there are many articles on performant Python specific to fractal generation around the web.
Another limitation, common to both the Python and OSL, is that the pattern is 2D . We’ll look into the 3D Mandelbulb later. For now, two simple tricks are to correlate the factor to an element’s scale on the z-axis
or to cast the 2D plane around a 3D form. In the case of a sphere, above right, the real and complex axis can be converted to spherical coordinates. A technique we’ll explore later on is to create a vertex group, weight the vertices in Python, then target that vertex group with modifiers.
Named after Gaston Julia, the Julia set differs from the Mandelbrot set insofar as the initial coordinate provided to the function is not set to a seed value the same way that
z was set to
To smooth the color gradient, we implement a technique posted on Stack Overflow by Mike Andrews and Per Alexandersson. The negative
ZN is supplied to the
exp function every iteration.
Fac records the sum.
The Wikipedia entry on the Julia set includes a gallery of interesting seeds. A handy method to find a seed is to use
rect to convert from polar to Cartesian coordinates. Set the radius to around
0.75, then cycle the angle
phi until a desired pattern emerges. This is also handy when animating a Julia set.
The workflow for the above in Python parallels that of the Mandelbrot set.
To make the dictionary returned by this function consistent with the above, we plug the original coordinate into the
'c' key, even though it’s referred to in the function as
z. This is so that the
visualize_fractal function can read the results of the
julia function below. The rest of the code to generate a Julia set matches that for the Mandelbrot set, so we’ll link to the gist here.
Pickover stalks, named after Clifford Pickover, introduce another way to visualize the algorithms we’ve already created above. Since the wispy strokes of this pattern require a high resolution, we won’t implement it in Python.
Instead of counting the number of iterations it takes for the
abs of a complex number to exceed the limit, we measure its absolute distance from an axis.
Between the horizontal and vertical distance, the minimum is always sought. After we’ve concluded the loop, this minimum is then divided by a number to govern the stroke weight of these stalks. The initial value of
TrapDist is an arbitrary large number.
Once we understand how to create Pickover stalks, we can incorporate them in to the same shaders as the original Julia and Mandelbrot sets.
The pattern lends itself well to the pinstriping look: when the
Fac generated by the above shader is greater than a lower-bound, it can be plugged into a
ColorRamp node in
RGB color mode with
Constant interpolation. A
Displacement node can be used to accentuate the strokes.
The Burning Ship
The Burning Ship fractal is straightforward, but is useful to illustrate the importance of controlling the input’s scale and offset. While the main burning ship figure is interesting, there are smaller, self-similar figures which include additional features, such as the cobwebby columns depicted below.
To zero in on the section of a fractal we want to see, we can use the vector
Mapping node or create our own. One option is to use use OSL’s
matrix data type. There’s more than one convention for the order in which these matrices are multiplied. Below, we use RST transformation order.
In OSL we can distinguish between
normals. These data types will not be affected by the
transform function in the same way. We cannot, however, distinguish between 2D texture coordinates and 3D coordinates. If we rotate a
UV coordinate, it will rotate from the bottom-left corner, not from the center,
Degrees can be converted to radians by multiplying them by π radians / 180 degrees (approximately
0.01745). The axis around which a point rotates should be normalized, but that is done within the shader, so we don’t have to worry about it when plugging in a value.
The rotation matrix above rotates clockwise (CW) when given a positive angle so as to match the behavior of OSL’s
Mapping node rotates
Texture coordinates counter-clockwise (CCW); the other types —
Normal — are rotated CW.
It may be preferable to not make a monolithic transform node at all. Translation is vector addition; scale, component-wise multiplication; rotation is covered by OSL’s
rotate. By letting each transformation stand on its own, we leave the transformation order (TRS, SRT, etc.) up to the composer.
Tricorn / Mandelbar
The tricorn or Mandelbar fractal raises the conjugate of z to a power.
The code above also demos how a Julia and Mandelbrot set can be toggled within the same shader using a Boolean
A 3D equivalent to the Mandelbrot set has long been the subject of experimentation. Daniel White and Paul Nylander have come up with a Mandelbulb by extending polar formula of Complex numbers to spherical formula. Left, Jonas Dichelle has a tutorial on how to create this with nodes in Blender’s new Eevee renderer. As with any 3D pattern, we can either emphasize a volume or the intersection of that volume with a surface.
A handy tool for the volumetric approach is the the relatively new Principled Volume shader node.
In order to recreate the Mandelbulb, we need to define spherical-Cartesian conversions. First, spherical to Cartesian,
clamp node is
min(max(value, lowerbound), upperbound). The above also groups
sin nodes together. Next, the reverse.
Beware: OSL’s global variables for angles are confusing.
M_PI_2 is π/2.
M_2PI is 2π.
M_2_PI is 2/π. In non-GPU programming languages, the Cartesian-spherical conversion should protect against a divide-by-zero exception.
Some versions of the formula use
cos of the inclination when converting to spherical and
acos when converting to Cartesian;
acos returns a value in [0, π]. Others use
sin of the inclination and
asin to reverse;
asin returns a value in [-π/2, π/2]. We keep with Barrallo’s “Expanding the Mandelbrot Set into Higher Dimensions,” and use
asin in the code to follow. For the spherical-Cartesian conversion, a precaution is to clamp the inclination and to ensure the radius is at least epsilon (ε), a small positive number.
Naming conventions may vary, with some formulae using theta (θ), phi (φ) and rho (ρ). This alternative adds uncertainty insofar as theta and phi may refer to azimuth and inclination or vice versa. Readers are encouraged to compare the world coordinates used by such formula against Blender’s.
To troubleshoot these conversions, we cut a corner out of a UV-sphere with the
Boolean modifier, then make a hue from the azimuth or the inclination. The sphere in the foreground on the left models inclination; in the background on the right, azimuth.
acos, a vector facing up on the z axis, (0, 0, 1), has an inclination of 0, and therefore is red. In addition, we’ve associated the radius with either value (background, right) or saturation (foreground, left).
The main formula entails raising the radius of
Z to an exponent; azimuth and inclination are multiplied by the exponent. Spherical coordinates are converted to Cartesian, then added to
Center. After this addition, the sum is used to update the spherical coordinates.
For Python, performance is even more of an issue than with 2D plots. The default grid count of 32 yields, worst case, 32, 768 cubes. Blender’s documentation offers a few best practices for writing scripts; particularly sound advice is to use list comprehensions and to time the script.
First, we spherical-Cartesian conversions in Python:
Then on to the main formula:
Just as we created a 2D grid for the Mandelbrot set, we create a 3D grid for the Mandelbulb:
As before, after these core functions are in place, the rest is a matter of composition. We won’t run through these, but the gist is here.
The same algorithm can be ported to a ‘Julia bulb’, but this simplistic approach doesn’t seem to have much traction compared to the approach that makes a 4D Julia set with quaternions. Inigo Quilez describes his own experiments with 3D Julia sets here.
Mixing A Fractal With Other Patterns
To demo but one approach for mixing a fractal with other patterns, we create a custom noise node; unlike the Cycles
Noise texture, this generates a
vector as an output. The greater aesthetic motive behind this approach is to reconcile the smoothness of noise with the exploration of roughness that motivated Mandelbrot in developing fractals.
The return value of the
noise function can be either a
vector or a
color. If we supply a
point and an additional float,
T, to the function, we can harness 4D noise. Simplex noise returns a value in the range [-1, 1] so we shift the
Fac to [0, 1]. The documentation includes a list of other noise alternatives:
"usimplex". The ‘u’ prefix stands for unsigned.
The output vector of a noise node can then contribute to the
Center input of a fractal shader. Noise in 4D makes it easy to animate, since all we need to do is key-frame the input
T. Alternatively, we could animate a point’s transformation, then then feed it into 3D noise.
For Python scripts that generate meshes, we can integrate them with Blender modifiers by setting vertex weights according to the
'fac' returned by our fractal functions.
First a new vertex group is created. Then the
visualize_fractal function can be updated. Then, a list of indices is added to the vertex group. The trick in this case is to ensure a fractal’s
'fac' can easily be tied back to the eight vertices of a cube; we simplify the function above by only creating a vertex group if no cubes have been pruned. With cube pruning, we don’t know how many vertices are in a mesh’s data; without, the list of vertices is of a predictable length.
If we don’t mind mutating the
results of the fractal in the process of visualizing it, we can update our
visualize_fractal function to remove data from the results list when pruning. This way, our abstract data corresponds with the visual data. In this case, the color of vertices is associated with the loops in a mesh’s data.
Unlike other programming languages, in Python,
pop — associated with the behavior of a
stack data structure — accepts an index as an argument. We could also use
remove; see this Stack Overflow thread for more. Since we are shifting elements, we loop through the list backwards.
For each cube in the fractal, there will be 24 loops, so that is the stride we’ll be working with. To create a color we could either use Python’s colorsys or Blender’s Color in
mathutils. Vertex weights are typically displayed with a color ramp that ranges from blue to red. To match this, we subtract
fac from one, then multiply by two-thirds. Not in the code above, we can also convert
z (or any other data returned by all fractal patterns) from a vector to a color.
If a material is supplied to the function, we can automate the process of adding an
Attribute input node and separating the color into HSV.
Instead of directly rendering vertex colors, we can convert a color back into a factor, [0, 1]. In the image above, such as factor is supplied to a noise texture, to a color ramp and to a displacement node.
If nothing else, learning complex numbers in 2D makes it considerably easier to understand how quaternions — which generalize complex numbers into 4D — make 3D rotations possible. The video from 3Blue1Brown, left, illustrates this connection.
A second take away for this article is how a generative algorithm is iterated. We began with a specific algorithm,
zn:= sq(z) + c, then studied variations. Each variation can be subdivided into how to generate new data, and then how to visualize it. The re-creation of these patterns is a starting point, not necessarily an end in itself. By studying past variations, we gain a sense of where opportunities for new visuals are to be found.
For more fractals, the Fractal Forums are a font of information. Paul Bourke’s website contains a gallery of variants: the bedouin, binocular, cactus, crown, lemon, log spiral, magnetic, secant sea and Zubieta fractal. Paul Nylander’s website does as well, particularly the page, “Hypercomplex Fractals.”
Beyond that, there are a number of fractal patterns which require no understanding of Complex numbers whatsoever, but emphasize the recursive, or iterative, process: Menger sponges, Koch curves, fractal flames, and so on. Noteworthy among 3D fractals is the Mandelbox, as it iterates spatial folding algorithms, specifically box- and sphere-folding. Fractal Lab demonstrates how the Mandelbox can be visualized in GLSL. This fractal has been explored by tree3d and DERBENDER at the Blender Artists Forums.