437 lines
13 KiB
Markdown
437 lines
13 KiB
Markdown
# Fractal Rendering Skill
|
|
|
|
## Use Cases
|
|
- Rendering self-similar mathematical structures: Mandelbrot/Julia sets (2D), Mandelbulb (3D), IFS fractals (Menger/Apollonian)
|
|
- Procedural textures or backgrounds requiring infinite detail
|
|
- Real-time generation of complex geometric visual effects (music visualization, sci-fi scenes, abstract art)
|
|
- Suitable for ShaderToy, demo scene, procedural content generation
|
|
|
|
## Core Principles
|
|
|
|
Fractal rendering is essentially **visualization of iterative systems**, falling into three categories:
|
|
|
|
### 1. Escape-Time Algorithm
|
|
Iterate `Z <- Z^2 + c`, count escape steps. Distance estimation by simultaneously tracking the derivative `Z'`:
|
|
```
|
|
Z <- Z^2 + c
|
|
Z' <- 2*Z*Z' + 1
|
|
d(c) = |Z|*log|Z| / |Z'|
|
|
```
|
|
|
|
### 2. Iterated Function System (IFS / KIFS)
|
|
Fold-sort-scale-offset iteration produces self-similar structures:
|
|
```
|
|
p = abs(p) // fold
|
|
sort p.xyz descending // sort
|
|
p = Scale * p - Offset * (Scale-1) // scale and offset
|
|
```
|
|
|
|
### 3. Spherical Inversion Fractals
|
|
`fract()` space folding + spherical inversion `p *= s/dot(p,p)`:
|
|
```
|
|
p = -1.0 + 2.0 * fract(0.5*p + 0.5)
|
|
k = s / dot(p, p)
|
|
p *= k; scale *= k
|
|
```
|
|
|
|
All 3D fractals are rendered via **Sphere Tracing (Ray Marching)**.
|
|
|
|
## Implementation Steps
|
|
|
|
### Step 1: Coordinate Normalization
|
|
```glsl
|
|
vec2 p = (2.0 * fragCoord - iResolution.xy) / iResolution.y;
|
|
```
|
|
|
|
### Step 2: 2D Mandelbrot Escape-Time Iteration
|
|
```glsl
|
|
float distanceToMandelbrot(in vec2 c) {
|
|
vec2 z = vec2(0.0);
|
|
vec2 dz = vec2(0.0);
|
|
float m2 = 0.0;
|
|
|
|
for (int i = 0; i < MAX_ITER; i++) {
|
|
if (m2 > BAILOUT * BAILOUT) break;
|
|
// Z' -> 2*Z*Z' + 1
|
|
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
|
|
z = vec2(z.x*z.x - z.y*z.y, 2.0*z.x*z.y) + c;
|
|
m2 = dot(z, z);
|
|
}
|
|
return 0.5 * sqrt(dot(z,z) / dot(dz,dz)) * log(dot(z,z));
|
|
}
|
|
```
|
|
|
|
### Step 3: Mandelbulb Distance Field (Spherical Coordinate Power-N)
|
|
```glsl
|
|
float mandelbulb(vec3 p) {
|
|
vec3 z = p;
|
|
float dr = 1.0;
|
|
float r;
|
|
|
|
for (int i = 0; i < FRACTAL_ITER; i++) {
|
|
r = length(z);
|
|
if (r > BAILOUT) 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;
|
|
}
|
|
```
|
|
|
|
### Step 4: Menger Sponge Distance Field (KIFS)
|
|
```glsl
|
|
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;
|
|
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;
|
|
}
|
|
```
|
|
|
|
### Step 5: Apollonian Distance Field (Spherical Inversion)
|
|
```glsl
|
|
vec4 orb; // 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);
|
|
float r2 = dot(p, p);
|
|
orb = min(orb, vec4(abs(p), r2));
|
|
float k = s / r2;
|
|
p *= k;
|
|
scale *= k;
|
|
}
|
|
return 0.25 * abs(p.y) / scale;
|
|
}
|
|
```
|
|
|
|
### Step 6: Ray Marching
|
|
```glsl
|
|
float rayMarch(vec3 ro, vec3 rd) {
|
|
float t = 0.01;
|
|
for (int i = 0; i < MAX_STEPS; i++) {
|
|
float precis = PRECISION * t;
|
|
float h = map(ro + rd * t);
|
|
if (h < precis || t > MAX_DIST) break;
|
|
t += h * FUDGE_FACTOR;
|
|
}
|
|
return (t > MAX_DIST) ? -1.0 : t;
|
|
}
|
|
```
|
|
|
|
### Step 7: Normal Calculation
|
|
```glsl
|
|
// 4-tap tetrahedral method (recommended)
|
|
vec3 calcNormal(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 & Lighting
|
|
```glsl
|
|
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);
|
|
|
|
vec3 brdf = vec3(0.4) * amb * ao + vec3(1.0) * diff * ao;
|
|
|
|
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
|
|
```glsl
|
|
void setupCamera(vec2 uv, vec3 ro, vec3 ta, float cr, out vec3 rd) {
|
|
vec3 cw = normalize(ta - ro);
|
|
vec3 cp = vec3(sin(cr), cos(cr), 0.0);
|
|
vec3 cu = normalize(cross(cw, cp));
|
|
vec3 cv = normalize(cross(cu, cw));
|
|
rd = normalize(uv.x * cu + uv.y * cv + 2.0 * cw);
|
|
}
|
|
```
|
|
|
|
## Complete Code Template
|
|
|
|
3D Apollonian fractal (spherical inversion type) with full ray marching pipeline, orbit trap coloring, and AO. Ready to run in ShaderToy.
|
|
|
|
```glsl
|
|
// Fractal Rendering — Apollonian (Spherical Inversion) Template
|
|
|
|
#define MAX_STEPS 200
|
|
#define MAX_DIST 30.0
|
|
#define PRECISION 0.001
|
|
#define INVERSION_ITER 8 // Tunable: 5-12
|
|
#define AA 1 // Tunable: 1=no AA, 2=4xSSAA
|
|
|
|
vec4 orb;
|
|
|
|
float map(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);
|
|
float r2 = dot(p, p);
|
|
orb = min(orb, vec4(abs(p), r2));
|
|
float k = s / r2;
|
|
p *= k;
|
|
scale *= k;
|
|
}
|
|
return 0.25 * abs(p.y) / scale;
|
|
}
|
|
|
|
float trace(vec3 ro, vec3 rd, float s) {
|
|
float t = 0.01;
|
|
for (int i = 0; i < MAX_STEPS; i++) {
|
|
float precis = PRECISION * t;
|
|
float h = map(ro + rd * t, s);
|
|
if (h < precis || t > MAX_DIST) break;
|
|
t += h;
|
|
}
|
|
return (t > MAX_DIST) ? -1.0 : t;
|
|
}
|
|
|
|
vec3 calcNormal(vec3 pos, float t, float s) {
|
|
float precis = PRECISION * t;
|
|
vec2 e = vec2(1.0, -1.0) * precis;
|
|
return normalize(
|
|
e.xyy * map(pos + e.xyy, s) +
|
|
e.yyx * map(pos + e.yyx, s) +
|
|
e.yxy * map(pos + e.yxy, s) +
|
|
e.xxx * map(pos + e.xxx, s));
|
|
}
|
|
|
|
vec3 render(vec3 ro, vec3 rd, float anim) {
|
|
vec3 col = vec3(0.0);
|
|
float t = trace(ro, rd, anim);
|
|
|
|
if (t > 0.0) {
|
|
vec4 tra = orb;
|
|
vec3 pos = ro + t * rd;
|
|
vec3 nor = calcNormal(pos, t, anim);
|
|
|
|
vec3 light1 = normalize(vec3(0.577, 0.577, -0.577));
|
|
vec3 light2 = normalize(vec3(-0.707, 0.0, 0.707));
|
|
float key = clamp(dot(light1, nor), 0.0, 1.0);
|
|
float bac = clamp(0.2 + 0.8 * dot(light2, nor), 0.0, 1.0);
|
|
float amb = 0.7 + 0.3 * nor.y;
|
|
float ao = pow(clamp(tra.w * 2.0, 0.0, 1.0), 1.2);
|
|
|
|
vec3 brdf = vec3(0.40) * amb * ao
|
|
+ vec3(1.00) * key * ao
|
|
+ vec3(0.40) * bac * ao;
|
|
|
|
vec3 rgb = vec3(1.0);
|
|
rgb = mix(rgb, vec3(1.0, 0.80, 0.2), clamp(6.0 * tra.y, 0.0, 1.0));
|
|
rgb = mix(rgb, vec3(1.0, 0.55, 0.0), pow(clamp(1.0 - 2.0*tra.z, 0.0, 1.0), 8.0));
|
|
|
|
col = rgb * brdf * exp(-0.2 * t);
|
|
}
|
|
return sqrt(col);
|
|
}
|
|
|
|
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
|
|
float time = iTime * 0.25;
|
|
float anim = 1.1 + 0.5 * smoothstep(-0.3, 0.3, cos(0.1 * iTime));
|
|
|
|
vec3 tot = vec3(0.0);
|
|
|
|
#if AA > 1
|
|
for (int jj = 0; jj < AA; jj++)
|
|
for (int ii = 0; ii < AA; ii++)
|
|
#else
|
|
int ii = 1, jj = 1;
|
|
#endif
|
|
{
|
|
vec2 q = fragCoord.xy + vec2(float(ii), float(jj)) / float(AA);
|
|
vec2 p = (2.0 * q - iResolution.xy) / iResolution.y;
|
|
|
|
vec3 ro = vec3(2.8*cos(0.1 + 0.33*time),
|
|
0.4 + 0.3*cos(0.37*time),
|
|
2.8*cos(0.5 + 0.35*time));
|
|
vec3 ta = vec3(1.9*cos(1.2 + 0.41*time),
|
|
0.4 + 0.1*cos(0.27*time),
|
|
1.9*cos(2.0 + 0.38*time));
|
|
float roll = 0.2 * cos(0.1 * time);
|
|
|
|
vec3 cw = normalize(ta - ro);
|
|
vec3 cp = vec3(sin(roll), cos(roll), 0.0);
|
|
vec3 cu = normalize(cross(cw, cp));
|
|
vec3 cv = normalize(cross(cu, cw));
|
|
vec3 rd = normalize(p.x*cu + p.y*cv + 2.0*cw);
|
|
|
|
tot += render(ro, rd, anim);
|
|
}
|
|
|
|
tot /= float(AA * AA);
|
|
fragColor = vec4(tot, 1.0);
|
|
}
|
|
```
|
|
|
|
## Common Variants
|
|
|
|
### 1. 2D Mandelbrot (Distance Estimation Coloring)
|
|
Pure 2D, no ray marching needed. Complex iteration + distance coloring.
|
|
```glsl
|
|
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
|
|
vec2 p = (2.0*fragCoord - iResolution.xy) / iResolution.y;
|
|
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
|
|
|
|
vec2 z = vec2(0.0), dz = vec2(0.0);
|
|
for (int i = 0; i < 300; i++) {
|
|
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);
|
|
fragColor = vec4(vec3(d), 1.0);
|
|
}
|
|
```
|
|
|
|
### 2. Mandelbulb Power-N
|
|
Spherical coordinate trigonometric functions; `POWER` parameter controls morphology.
|
|
```glsl
|
|
#define POWER 8.0 // Tunable: 2-16
|
|
#define FRACTAL_ITER 4 // Tunable: 2-8
|
|
|
|
float mandelbulbDE(vec3 p) {
|
|
vec3 z = p;
|
|
float dr = 1.0, 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)
|
|
`abs()` folding + conditional sorting, regular geometric fractal.
|
|
```glsl
|
|
#define SCALE 3.0
|
|
#define OFFSET vec3(0.92858,0.92858,0.32858)
|
|
#define IFS_ITER 7
|
|
|
|
float mengerDE(vec3 z) {
|
|
z = abs(1.0 - mod(z, 2.0));
|
|
float d = 1000.0;
|
|
for (int n = 0; n < IFS_ITER; n++) {
|
|
z = abs(z);
|
|
if (z.x < z.y) z.xy = z.yx;
|
|
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
|
|
Quaternion `Z <- Z^2 + c` (4D), with fixed `c` parameter; visualized by taking a 3D slice.
|
|
```glsl
|
|
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, mz2 = dot(z, z);
|
|
for (int i = 0; i < 11; i++) {
|
|
md2 *= 4.0 * mz2;
|
|
z = qsqr(z) + c;
|
|
mz2 = dot(z, z);
|
|
if (mz2 > 4.0) break;
|
|
}
|
|
return 0.25 * sqrt(mz2 / md2) * log(mz2);
|
|
}
|
|
// Animated 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)
|
|
`abs(p)/dot(p,p) + offset` iteration, weighted accumulation produces a density field.
|
|
```glsl
|
|
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) {
|
|
float mag = dot(p, p);
|
|
p = abs(p) / mag + vec3(-0.5, -0.4, -1.5); // Tunable: offset values
|
|
float w = exp(-float(i) / 7.0);
|
|
accum += w * exp(-strength * pow(abs(mag - prev), 2.3));
|
|
tw += w;
|
|
prev = mag;
|
|
}
|
|
return max(0.0, 5.0 * accum / tw - 0.7);
|
|
}
|
|
```
|
|
|
|
## Performance & Composition
|
|
|
|
### Performance Tips
|
|
- Core bottleneck: outer ray marching x inner fractal iteration (e.g., `200 x 8 = 1600` map calls per pixel)
|
|
- Reduce `MAX_STEPS` to 60-100, compensate with fudge factor 0.7-0.9
|
|
- Hit threshold `precis = 0.001 * t` relaxes with distance
|
|
- Fractal iteration: break immediately when `|z|^2 > bailout`
|
|
- Reducing iterations from 8 to 4-5 has minimal visual impact
|
|
- Use 4-tap normals instead of 6-tap to save 33%
|
|
- Use AA=1 during development, AA=2 for release (AA=3 = 9x overhead)
|
|
- Avoid `pow()` inside loops; manually expand for low powers
|
|
|
|
### Composition Techniques
|
|
- **Volumetric light**: accumulate `exp(-10.0 * h)` during ray march for god rays
|
|
- **Tone Mapping**: ACES + sRGB gamma for handling high-frequency detail
|
|
- **Transparent refraction**: negative distance field reverse ray march + Beer's law absorption
|
|
- **Orbit Trap coloring**: map trap values to HSV or emissive colors
|
|
- **Soft shadows**: ray march toward light, accumulate `min(k * h / t)` for soft shadows
|
|
|
|
## Further Reading
|
|
|
|
For complete step-by-step tutorials, mathematical derivations, and advanced usage, see [reference](../reference/fractal-rendering.md)
|