18 KiB
Fractal Rendering — Detailed Reference
This document is a detailed supplement to SKILL.md, containing prerequisites, step-by-step explanations, mathematical derivations, variant descriptions, in-depth performance analysis, and complete combination example code.
Prerequisites
- GLSL Basics: uniform, varying, built-in functions (
dot,length,normalize,abs,fract) - Complex Number Arithmetic: representing complex numbers as
vec2, multiplication(a+bi)(c+di) = (ac-bd, ad+bc) - Vector Math: dot product, cross product, matrix transforms
- Ray Marching Basics (required for 3D fractals): stepping along a ray, using distance fields for collision detection
- Coordinate Normalization: mapping pixel coordinates to the
[-1, 1]range
Core Principles in Detail
The essence of fractal rendering is visualization of iterative systems. Core algorithm patterns fall into three categories:
1. Escape-Time Algorithm
For each point c on the complex plane, repeatedly iterate Z <- Z^2 + c, counting the number of steps needed for Z to escape (|Z| > R). More steps means closer to the fractal boundary.
Distance Estimation computes the precise distance from a point to the fractal by simultaneously tracking the derivative Z':
Z <- Z^2 + c (value iteration)
Z' <- 2*Z*Z' + 1 (derivative iteration)
d(c) = |Z|*log|Z| / |Z'| (Hubbard-Douady potential function)
Distance estimation produces smoother coloring than pure escape-time step counting, and is a prerequisite for ray marching in 3D fractals.
2. Iterated Function System (IFS)
Apply a set of transforms (folding abs(), scaling Scale, offset Offset) to points in space, iterating repeatedly to produce self-similar structures. Core steps of KIFS (Kaleidoscopic IFS) commonly used in 3D:
p = abs(p) // Fold (symmetrize)
sort p.xyz descending // Sort (select symmetry axis)
p = Scale * p - Offset * (Scale-1) // Scale and offset
3. Spherical Inversion Fractal
Apollonian-type fractals use fract() for space folding + spherical inversion p *= s/dot(p,p):
p = -1.0 + 2.0 * fract(0.5*p + 0.5) // Fold space to [-1,1]
r^2 = dot(p, p)
k = s / r^2 // Inversion factor
p *= k; scale *= k // Spherical inversion
All 3D fractals are rendered using Sphere Tracing (Ray Marching): stepping along the view ray by the distance field value at each step, until close enough to the surface.
Implementation Steps in Detail
Step 1: Coordinate Normalization
What: Map pixel coordinates to standard coordinates centered on the screen with aspect ratio correction.
Why: All fractal calculations must be performed in mathematical space, independent of pixel resolution.
vec2 p = (2.0 * fragCoord - iResolution.xy) / iResolution.y;
// p now has y range [-1,1], x scaled by aspect ratio
Step 2: 2D Fractal — Mandelbrot Escape-Time Iteration
What: For each pixel point as complex number c, iterate Z <- Z^2 + c while tracking the derivative.
Why: Escape time produces fractal structure; derivative tracking enables distance estimation coloring.
float distanceToMandelbrot(in vec2 c) {
vec2 z = vec2(0.0);
vec2 dz = vec2(0.0); // Derivative
float m2 = 0.0;
for (int i = 0; i < MAX_ITER; i++) {
if (m2 > BAILOUT * BAILOUT) break;
// Z' -> 2*Z*Z' + 1 (complex derivative chain rule)
dz = 2.0 * vec2(z.x*dz.x - z.y*dz.y,
z.x*dz.y + z.y*dz.x) + vec2(1.0, 0.0);
// Z -> Z^2 + c (complex squaring)
z = vec2(z.x*z.x - z.y*z.y, 2.0*z.x*z.y) + c;
m2 = dot(z, z);
}
// Distance estimation: d(c) = |Z|*log|Z| / |Z'|
return 0.5 * sqrt(dot(z,z) / dot(dz,dz)) * log(dot(z,z));
}
Step 3: 3D Fractal — Distance Field Function (Mandelbulb Example)
What: Implement the Mandelbulb power-N iteration using spherical coordinates, returning a distance estimate.
Why: 3D fractals cannot be directly colored via escape-time on pixels; they require distance fields for ray marching.
float mandelbulb(vec3 p) {
vec3 z = p;
float dr = 1.0; // Derivative (distance scaling factor)
float r;
for (int i = 0; i < FRACTAL_ITER; i++) {
r = length(z);
if (r > BAILOUT) break;
// Convert to spherical coordinates
float theta = atan(z.y, z.x);
float phi = asin(z.z / r);
// Derivative: dr -> power * r^(power-1) * dr + 1
dr = pow(r, POWER - 1.0) * dr * POWER + 1.0;
// z -> z^power + p (spherical coordinate exponentiation)
r = pow(r, POWER);
theta *= POWER;
phi *= POWER;
z = r * vec3(cos(theta)*cos(phi),
sin(theta)*cos(phi),
sin(phi)) + p;
}
// Distance estimation
return 0.5 * log(r) * r / dr;
}
Step 4: 3D Fractal — IFS Distance Field (Menger Sponge Example)
What: Construct a KIFS fractal distance field through fold-sort-scale-offset iteration.
Why: IFS fractals produce self-similar structures through spatial transforms rather than numerical iteration; distance is tracked via Scale^(-n) scaling.
float mengerDE(vec3 z) {
z = abs(1.0 - mod(z, 2.0)); // Infinite tiling
float d = 1000.0;
for (int n = 0; n < IFS_ITER; n++) {
z = abs(z); // Fold
if (z.x < z.y) z.xy = z.yx; // Sort
if (z.x < z.z) z.xz = z.zx;
if (z.y < z.z) z.yz = z.zy;
z = SCALE * z - OFFSET * (SCALE - 1.0); // Scale + offset
if (z.z < -0.5 * OFFSET.z * (SCALE - 1.0))
z.z += OFFSET.z * (SCALE - 1.0);
d = min(d, length(z) * pow(SCALE, float(-n) - 1.0));
}
return d - 0.001;
}
Step 5: 3D Fractal — Spherical Inversion Distance Field (Apollonian Type)
What: Construct an Apollonian fractal using fract folding + spherical inversion iteration, while recording orbit traps.
Why: Spherical inversion p *= s/dot(p,p) produces sphere packing structures; orbit traps provide color and AO information.
vec4 orb; // Global orbit trap
float apollonianDE(vec3 p, float s) {
float scale = 1.0;
orb = vec4(1000.0);
for (int i = 0; i < INVERSION_ITER; i++) {
p = -1.0 + 2.0 * fract(0.5 * p + 0.5); // Fold space to [-1,1]
float r2 = dot(p, p);
orb = min(orb, vec4(abs(p), r2)); // Record orbit trap
float k = s / r2; // Inversion factor
p *= k;
scale *= k;
}
return 0.25 * abs(p.y) / scale;
}
Step 6: Ray Marching (Sphere Tracing)
What: Step along the ray direction, advancing by the distance field value at each step, until hitting the surface.
Why: The distance field guarantees safe stepping (won't pass through the surface), and is the standard method for rendering implicit 3D fractals.
float rayMarch(vec3 ro, vec3 rd) {
float t = 0.01;
for (int i = 0; i < MAX_STEPS; i++) {
float precis = PRECISION * t; // Relax precision with distance
float h = map(ro + rd * t);
if (h < precis || t > MAX_DIST) break;
t += h * FUDGE_FACTOR; // fudge < 1.0 improves safety
}
return (t > MAX_DIST) ? -1.0 : t;
}
Step 7: Normal Calculation (Finite Differences)
What: Sample the distance field gradient around the hit point as the surface normal.
Why: Implicit surfaces have no analytical normals and require numerical approximation. Tetrahedral sampling (4-tap) saves 1/3 of the cost compared to central differences (6-tap).
// 6-tap central difference method (more intuitive)
vec3 calcNormal_6tap(vec3 pos) {
vec2 e = vec2(0.001, 0.0);
return normalize(vec3(
map(pos + e.xyy) - map(pos - e.xyy),
map(pos + e.yxy) - map(pos - e.yxy),
map(pos + e.yyx) - map(pos - e.yyx)));
}
// 4-tap tetrahedral method (more efficient, recommended)
vec3 calcNormal_4tap(vec3 pos, float t) {
float precis = 0.001 * t;
vec2 e = vec2(1.0, -1.0) * precis;
return normalize(
e.xyy * map(pos + e.xyy) +
e.yyx * map(pos + e.yyx) +
e.yxy * map(pos + e.yxy) +
e.xxx * map(pos + e.xxx));
}
Step 8: Shading and Lighting
What: Compute Lambertian diffuse + ambient + AO for hit surfaces.
Why: Lighting gives 3D fractals depth and material quality. Orbit trap values (orb) can serve both as color mapping and as simple AO.
vec3 shade(vec3 pos, vec3 nor, vec3 rd, vec4 trap) {
vec3 light1 = normalize(LIGHT_DIR);
float diff = clamp(dot(light1, nor), 0.0, 1.0);
float amb = 0.7 + 0.3 * nor.y;
float ao = pow(clamp(trap.w * 2.0, 0.0, 1.0), 1.2); // Orbit trap AO
vec3 brdf = vec3(0.4) * amb * ao // Ambient
+ vec3(1.0) * diff * ao; // Diffuse
// Map material color from orbit trap
vec3 rgb = vec3(1.0);
rgb = mix(rgb, vec3(1.0, 0.8, 0.2), clamp(6.0*trap.y, 0.0, 1.0));
rgb = mix(rgb, vec3(1.0, 0.55, 0.0), pow(clamp(1.0-2.0*trap.z, 0.0, 1.0), 8.0));
return rgb * brdf;
}
Step 9: Camera Setup
What: Build a look-at camera matrix, converting pixel coordinates to 3D ray directions.
Why: All 3D fractal ray marching requires a unified camera framework to generate rays.
void setupCamera(vec2 uv, vec3 ro, vec3 ta, float cr,
out vec3 rd) {
vec3 cw = normalize(ta - ro); // forward
vec3 cp = vec3(sin(cr), cos(cr), 0.0); // roll
vec3 cu = normalize(cross(cw, cp)); // right
vec3 cv = normalize(cross(cu, cw)); // up
rd = normalize(uv.x * cu + uv.y * cv + 2.0 * cw); // FOV ~ 2.0
}
Common Variants in Detail
1. 2D Mandelbrot (Distance Estimation Coloring)
Difference from base version (3D Apollonian): pure 2D computation, no ray marching needed, uses complex iteration + distance coloring.
// Replace entire mainImage
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 p = (2.0*fragCoord - iResolution.xy) / iResolution.y;
// Animated zoom
float tz = 0.5 - 0.5*cos(0.225*iTime);
float zoo = pow(0.5, 13.0*tz);
vec2 c = vec2(-0.05, 0.6805) + p * zoo; // Tunable: zoom center point
// Iteration
vec2 z = vec2(0.0), dz = vec2(0.0);
for (int i = 0; i < 300; i++) { // Tunable: iteration count
if (dot(z,z) > 1024.0) break;
dz = 2.0*vec2(z.x*dz.x-z.y*dz.y, z.x*dz.y+z.y*dz.x) + vec2(1.0,0.0);
z = vec2(z.x*z.x-z.y*z.y, 2.0*z.x*z.y) + c;
}
float d = 0.5*sqrt(dot(z,z)/dot(dz,dz))*log(dot(z,z));
d = clamp(pow(4.0*d/zoo, 0.2), 0.0, 1.0); // Tunable: 0.2 controls contrast
fragColor = vec4(vec3(d), 1.0);
}
2. Mandelbulb Power-N (3D Spherical Coordinate Fractal)
Difference from base version: uses spherical coordinate trigonometric functions instead of spherical inversion, with a tunable POWER parameter controlling the fractal shape.
#define POWER 8.0 // Tunable: 2-16, higher = more complex structure
#define FRACTAL_ITER 4 // Tunable: 2-8, more = more detail
float mandelbulbDE(vec3 p) {
vec3 z = p;
float dr = 1.0;
float r;
for (int i = 0; i < FRACTAL_ITER; i++) {
r = length(z);
if (r > 2.0) break;
float theta = atan(z.y, z.x);
float phi = asin(z.z / r);
dr = pow(r, POWER - 1.0) * dr * POWER + 1.0;
r = pow(r, POWER);
theta *= POWER;
phi *= POWER;
z = r * vec3(cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi)) + p;
}
return 0.5 * log(r) * r / dr;
}
3. Menger Sponge (KIFS Folding Type)
Difference from base version: uses abs() folding + conditional sorting instead of spherical inversion, producing regular geometric fractals.
#define SCALE 3.0 // Tunable: scaling factor, 2.0-4.0
#define OFFSET vec3(0.92858,0.92858,0.32858) // Tunable: offset vector, changes shape
#define IFS_ITER 7 // Tunable: iteration count
float mengerDE(vec3 z) {
z = abs(1.0 - mod(z, 2.0)); // Infinite tiling
float d = 1000.0;
for (int n = 0; n < IFS_ITER; n++) {
z = abs(z);
if (z.x < z.y) z.xy = z.yx; // Conditional sorting
if (z.x < z.z) z.xz = z.zx;
if (z.y < z.z) z.yz = z.zy;
z = SCALE * z - OFFSET * (SCALE - 1.0);
if (z.z < -0.5*OFFSET.z*(SCALE-1.0))
z.z += OFFSET.z*(SCALE-1.0);
d = min(d, length(z) * pow(SCALE, float(-n)-1.0));
}
return d - 0.001;
}
4. Quaternion Julia Set
Difference from base version: uses quaternion algebra Z <- Z^2 + c (4D), Julia sets use a fixed c parameter instead of per-point c, visualized by taking 3D cross-sections.
// Quaternion squaring
vec4 qsqr(vec4 a) {
return vec4(a.x*a.x - a.y*a.y - a.z*a.z - a.w*a.w,
2.0*a.x*a.y, 2.0*a.x*a.z, 2.0*a.x*a.w);
}
float juliaDE(vec3 p, vec4 c) {
vec4 z = vec4(p, 0.0);
float md2 = 1.0;
float mz2 = dot(z, z);
for (int i = 0; i < 11; i++) { // Tunable: iteration count
md2 *= 4.0 * mz2; // |dz| -> 2*|z|*|dz|
z = qsqr(z) + c; // z -> z^2 + c
mz2 = dot(z, z);
if (mz2 > 4.0) break;
}
return 0.25 * sqrt(mz2 / md2) * log(mz2);
}
// Animated Julia parameter c:
// vec4 c = 0.45*cos(vec4(0.5,3.9,1.4,1.1) + time*vec4(1.2,1.7,1.3,2.5)) - vec4(0.3,0,0,0);
5. Minimal IFS Field (2D, No Ray Marching)
Difference from base version: pure 2D implementation, only ~20 lines of code, using abs(p)/dot(p,p) + offset for iteration, producing a density field through weighted accumulation.
float field(vec3 p) {
float strength = 7.0 + 0.03 * log(1.e-6 + fract(sin(iTime) * 4373.11));
float accum = 0.0, prev = 0.0, tw = 0.0;
for (int i = 0; i < 32; ++i) { // Tunable: iteration count
float mag = dot(p, p);
p = abs(p) / mag + vec3(-0.5, -0.4, -1.5); // Tunable: offset values change shape
float w = exp(-float(i) / 7.0); // Tunable: 7.0 controls decay
accum += w * exp(-strength * pow(abs(mag - prev), 2.3));
tw += w;
prev = mag;
}
return max(0.0, 5.0 * accum / tw - 0.7);
}
// Sample field() directly on fragCoord as brightness/color
Performance Optimization Details
Bottleneck Analysis
The core bottleneck in fractal rendering is nested loops: outer ray marching steps x inner fractal iterations. A single pixel may execute 200 steps x 8 iterations = 1600 distance field evaluations.
Optimization Techniques
1. Reduce Ray Marching Steps
Lower MAX_STEPS from 200 to 60-100, compensating precision loss with a fudge factor (0.7-0.9).
t += h * 0.7; // Fudge factor < 1.0, allows larger steps but reduces penetration risk
2. Adaptive Precision
Relax the collision threshold as distance increases; far objects don't need pixel-level precision.
float precis = 0.001 * t; // Precision grows linearly with distance
3. Early Exit
In fractal iteration, break immediately once |z|^2 > bailout.
if (m2 > 4.0) break; // Don't continue useless iterations
4. Reduce Iteration Count
Fractal iteration counts (INVERSION_ITER, IFS_ITER) reduced from 8 to 4-5 have minimal visual impact but significant performance gains.
5. Use 4-Tap Instead of 6-Tap for Normals
The tetrahedral method requires only 4 map() calls instead of 6, saving 33% normal computation cost.
6. AA Downgrade
Use #define AA 1 during development, switch to AA 2 for release. AA 3 has massive performance impact (9x overhead).
7. Distance Field Scaling
For non-unit-sized fractals, scale the space first then scale the distance value to avoid precision issues.
float z1 = 2.0;
return mandelbulb(p / z1) * z1;
8. Avoid pow() Inside Loops
pow(r, power) in Mandelbulb is expensive; low powers (e.g., 2, 3) can be manually expanded instead.
Combination Suggestions
1. Fractal + Volumetric Lighting
Accumulate scattered light passing through fractal gaps during ray marching, producing "god rays" effects.
// Accumulate additionally in ray march loop
float glow = 0.0;
for (...) {
float h = map(ro + rd*t);
glow += exp(-10.0 * h); // Closer to surface = larger contribution
...
}
col += glowColor * glow * 0.01;
2. Fractal + Post-Processing (Tone Mapping / FXAA)
3D fractals have rich high-frequency detail, prone to aliasing. Use ACES Tone Mapping + sRGB correction + FXAA post-processing.
// ACES tone mapping
vec3 aces_approx(vec3 v) {
v = max(v, 0.0) * 0.6;
float a=2.51, b=0.03, c=2.43, d=0.59, e=0.14;
return clamp((v*(a*v+b))/(v*(c*v+d)+e), 0.0, 1.0);
}
col = aces_approx(col);
col = pow(col, vec3(1.0/2.4)); // sRGB gamma
3. Fractal + Transparent Refraction (Multi-Bounce Refraction)
Used for "crystal ball" effects on volumetric fractals like Mandelbulb. Uses negative distance fields for reverse ray marching inside, combined with Beer's law absorption.
// Invert distance field for interior stepping
float dfactor = isInside ? -1.0 : 1.0;
float d = dfactor * map(ro + rd * t);
// Beer's law light absorption
ragg *= exp(-st * beer); // beer = negative color vector
// Refraction direction
vec3 refr = refract(rd, sn, isInside ? 1.0/ior : ior);
4. Fractal + Orbit Trap Texture Mapping
Orbit trap values can be mapped to HSV color space for rich coloring, or mapped as self-emission for glowing fractal effects.
vec3 hsv2rgb(vec3 c) {
vec4 K = vec4(1.0, 2.0/3.0, 1.0/3.0, 3.0);
vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}
// Map orbit trap to HSV
vec3 col = hsv2rgb(vec3(trap.x * 0.5, 0.9, 0.8));
5. Fractal + Soft Shadow
Perform an additional ray march from the fractal surface toward the light source, accumulating the minimum h/t ratio to generate soft shadows.
float softshadow(vec3 ro, vec3 rd, float mint, float k) {
float res = 1.0;
float t = mint;
for (int i = 0; i < 64; i++) {
float h = map(ro + rd*t);
res = min(res, k * h / t); // Larger k = harder shadows
if (res < 0.001) break;
t += clamp(h, 0.01, 0.5);
}
return clamp(res, 0.0, 1.0);
}