Visualizing Complex Numbers in Blender
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.
Background
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 0.375.
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
atan2,ceil,floor,fractandsqrt. 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 astructas 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 thecmathlibrary.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 vectorsaandb,a @ bis the dot product. For matrixmand vectorv,m @ vapplies the affine transformationmtov. For quaternionqand vectorv,q @ vrotatesvbyq.
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 abssq, abs, phase, conjugate and reciprocal.
Abs Sq / Modulus Sq / R Sq
Depends on: add, multiply, separate xyzNeeded for: abs, reciprocal
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, a.dot(a) or Vector.dot(a, a). This can be used to avoid the sqrt implicit in magnitude.
Abs / Modulus / R
Depends on: abssq, sqrt (or power)
Needed for: mandelbrot (and variants), polar
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, x and 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 z is 0.0.
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
Depends on: atan2 (or atan), separate xyzNeeded for: polar
atan2 takes the y or imag component as its first argument; the x or 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 trunc, not 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 +/-Infinity or NaN (JavaScript). atan(0.0)returns 0.0.
Since the 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
With Booleans, 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 phase above: float a = mix(0.0, M_PI, real < 0.0); can be simplified to float a = M_PI * (real < 0.0);.
Conjugate
Depends on: combine xyz, multiply, separate xyz
conjugate is usually too trivial to implement. It can, however, be an extra layer of abstraction to clarify reciprocal below.
In math notation, the conjugate of z has a horizontal bar over it. The conjugate can also be used to find the absolute,
as will become clearer after we address complex multiplication.
Reciprocal / Inverse
Depends on: abssqNeeded for: div
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.
Conversions
In this section we cover the conversions polar and rect.
Polar
Depends on: abs, phaseNeeded for: log
Rect
Depends on: combine xyz, cosine, sine, multiplyNeeded for: exp
Binary Operations, Pt. 1
In this section, we cover add, sub, mult and div. We also define pow, though it depends on exp and log, which have not yet been defined.
Add
Depends on: add, combine xyz, separate xyz
Needed for: mandelbrot (and variants), mobius
Sub
Depends on: combine xyz, subtract, separate xyz
Mult
Depends on: add, combine xyz, multiply, separate xyzNeeded for: div, mobius, pow
Div
Depends on: mult, reciprocalNeeded for: mobius
Pow
Depends on: exp, log, multNeeded for: mandelbrot (and variants)
Unary Operations, Pt. 2
This section covers exp, log, cos and sin.
Exp
Depends on: power (exp), rect, separate xyzNeeded for: pow
Log
Depends on: logarithm (ln), polarNeeded for: pow
Cos
Depends on: combine xyz, cos, cosh, multiply, separate xyz, sin, sinh
To implement the above in Cycles nodes, we need hyperbolic cosine (cosh)
and hyperbolic sine (sinh).
To offer a sense of the pattern created by cos, we use a diagnostic texture on a plane to animate the following:
Sin
Depends on: combine xyz, cos, cosh, multiply, separate xyz, sin, sinh
2D Patterns
Once we start with fractals, we’ll stop using Cycles nodes. This is because fractal algorithms require us to iterate a core algorithm.
Mobius Transformation
The Mobius transformation is named after August Möbius, the same man behind the Mobius strip, and is key to math behind projections.
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.
One way we can visualize this in Python is to create a lattice; we can then add a lattice modifier to a mesh object. First, we import some tools.
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.
Depending on which rotation mode is preferred, the quaternion can be converted using either to_axis_angle or to_euler. We include options for all three rotation modes above, but only one is needed.
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.
Mandelbrot Set
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 (2, 3, 4…), e is a positive real number (2.0, 2.1, 2.2…), e is a signed real number (-2.0, -1.5, 0.0, 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.
OSL code is below, but a demo of a node-based approach can be found in the video “Fractals in Blender (Cycles, Volumetrics)” from Kostack Studio.
We use floats in place of Booleans; when ExcludeUpper is equal to zero, it is false. When 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 complex.
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 for loop.
The list 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 offset and 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:
A 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 mandelbrot function:
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 float_lerp.
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.
Julia Set
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 center previously.
To smooth the color gradient, we implement a technique posted on Stack Overflow by Mike Andrews and Per Alexandersson. The negative Abs of 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
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 points, vectors and 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, (0.5, 0.5).
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 rotate. The Mapping node rotates Texture coordinates counter-clockwise (CCW); the other types — Point, Vector and 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 float.
3D Patterns
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.
Spherical-Cartesian Conversion
In order to recreate the Mandelbulb, we need to define spherical-Cartesian conversions. First, spherical to Cartesian,
The clamp node is min(max(value, lowerbound), upperbound). The above also groups cos and 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.
Using 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 Mandelbulb
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
Noise
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 float, 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: "cell", "gabor", "hash", "perlin", "uperlin", "usimplex". The ‘u’ prefix stands for unsigned.
The output vector of a noise node can then contribute to the Z or 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.
Other combinations include combining the output Facs of two fractals with logic gates (such as AND, OR, XOR) or supplying the final iteration of ZN as an input coordinate to a texture.
Vertex Weighting
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.
We can then specify the vertex group in any number of modifiers which support the option. In the animation above, a wave modifier is used.
Vertex Colors
Another technique is to paint vertex colors. We can then pass these colors into a material via an Attribute node.
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 del or 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.
Conclusion
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.
