レイマーチングで使える複雑めな距離関数を収集していきます。
球や直方体など単純な図形に関しては以下にまとまっているので、ここでは触れません。
Inigo Quilez :: fractals, computer graphics, mathematics, shaders, demoscene and more
Recursive Tetrahedron
http://glslsandbox.com/e#55834.0
四面体のIFSによる繰り返しです。最終的に四面体の頂点に対して距離を計算しています。
#define ITERATIONS 8
float deRecursiveTetrahedron(vec3 p, vec3 offset, float scale) {
vec4 z = vec4(p, 1.0);
for (int i = 0; i < ITERATIONS; i++) {
if (z.x + z.y < 0.0) z.xy = -z.yx;
if (z.x + z.z < 0.0) z.xz = -z.zx;
if (z.y + z.z < 0.0) z.zy = -z.yz;
z *= scale;
z.xyz -= offset * (scale - 1.0);
}
return (length(z.xyz) - 1.5) / z.w;
}
float de(vec3 p) {
return deRecursiveTetrahedron(p, vec3(1.0), 2.0);
}
参考
Menger Sponge
メンガーのスポンジです。立方体から再帰的にクロス状の図形との論理差を取ることで穴を開けていきます。
http://glslsandbox.com/e#55835.0
float sdCross(vec3 p, float c) {
p = abs(p);
float dxy = max(p.x, p.y);
float dyz = max(p.y, p.z);
float dxz = max(p.x, p.z);
return min(dxy, min(dyz, dxz)) - c;
}
float sdBox(vec3 p, vec3 b) {
p = abs(p) - b;
return length(max(p, 0.0)) + min(max(p.x, max(p.y, p.z)), 0.0);
}
#define ITERATIONS 5
float deMengerSponge1(vec3 p, float scale, float width) {
float d = sdBox(p, vec3(1.0));
float s = 1.0;
for (int i = 0; i < ITERATIONS; i++) {
vec3 a = mod(p * s, 2.0) - 1.0;
s *= scale;
vec3 r = 1.0 - scale * abs(a);
float c = sdCross(r, width) / s;
d = max(d, c);
}
return d;
}
float de(vec3 p) {
return deMengerSponge1(p, 3.0, 1.0);
}
メンガーのスポンジの別の実装です。こちらはIFSを使ったシンプルな実装になっています。
http://glslsandbox.com/e#55836.0
// https://gam0022.net/blog/2019/06/25/unity-raymarching/
#define ITERATIONS 5
float deMengerSponge2(vec3 p, vec3 offset, float scale) {
vec4 z = vec4(p, 1.0);
for (int i = 0; i < ITERATIONS; i++) {
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.xyz -= offset * (scale - 1.0);
if (z.z < -0.5 * offset.z * (scale - 1.0))
z.z += offset.z * (scale - 1.0);
}
return (length(max(abs(z.xyz) - vec3(1.0, 1.0, 1.0), 0.0))) / z.w;
}
float de(vec3 p) {
return deMengerSponge2(p, vec3(1.0), 3.0);
}
参考
- GLSLでIFSシェルピンスキーのカーペット - Qiita
- Inigo Quilez :: fractals, computer graphics, mathematics, shaders, demoscene and more
- Unity×レイマーチングによる映像制作の実践手法 | gam0022.net
Hexagonal Tiling
http://glslsandbox.com/e#55837.0
六角形のタイリングです。六角形ではグリッド状にきれいに分割することができないので、単純にmodを用いた繰り返しを使用してもうまくいきません。そのため、六角形のタイリングを2つに分けてそれぞれ距離を計算し、論理和で結合しています。
float sdHex(vec2 p, float h) {
vec3 k = vec3(-0.8660254, 0.57735, 0.5);
p = abs(p);
p -= 2.0 * min(dot(k.xz, p), 0.0) * k.xz;
return length(p - vec2(clamp(p.x, -k.y * h, k.y * h), h)) * sign(p.y - h);
}
// SQRT3 = sqrt(3.0)
#define SQRT3 1.73205080757
// 0.0 <= scale <= 1.0
float deHexTiling(vec2 p, float radius, float scale) {
vec2 rep = vec2(2.0 * SQRT3, 2.0) * radius;
vec2 p1 = mod(p, rep) - rep * 0.5;
vec2 p2 = mod(p + 0.5 * rep, rep) - rep * 0.5;
return min(
sdHex(p1.xy, scale * radius),
sdHex(p2.xy, scale * radius)
);
}
float de(vec3 p) {
return max(deHexTiling(p.zx, 1.0, 0.9), p.y);
}
Quaternion Mandelbrot Set
http://glslsandbox.com/e#55838.0
マンデルブロ集合のクォータニオン版です。マンデルブロ集合は以下の漸化式で発散しない複素数(クォータニオン)$c$の集合になります。
\begin{align}
z_{0} &= 0\\
z_{n+1} &= z_{n}^{2} + c
\end{align}
この漸化式では$c$で発散するかどうかしかわかりません。レイマーチングでは物体表面までの距離が必要になりますが、Distance Estimationという手法により距離の下限値は以下の式で求まります。
DE = 0.5r(\ln{r}) / dr
ここで$r$はescape time length、$dr$はlength of running derivativeになり、以下により求まります。
(私も詳しく説明できないので、詳細に知りたい方はこの記事を読んでください)
f_{n}(c) = f_{n-1}^{2}(c) + c, f_{0} = 0\\
f'_{n}(c) = 2f_{n - 1}(c)f'_{n - 1}(c) + 1\\
r = |f_{n}(c)|\\
dr = |f'_{n}(c)|
これを愚直に実装したものが以下のようになります。漸化式の積がクォータニオンのものになっていることに注意してください。
vec4 qmul(vec4 a, vec4 b) {
return vec4(
a.x * b.x - a.y * b.y - a.z * b.z - a.w * b.w,
a.x * b.y + a.y * b.x - a.z * b.w + a.w * b.z,
a.x * b.z + a.y * b.w + a.z * b.x - a.w * b.y,
a.x * b.w - a.y * b.z + a.z * b.y + a.w * b.x
);
}
#define ITERATIONS 16
float deQuaternionMandelbrot(vec4 p) {
vec4 z = vec4(0.0);
vec4 dz = vec4(0.0);
vec4 pz, pdz;
float r, dr;
for (int i = 0; i < ITERATIONS; i++) {
pz = z;
z = qmul(pz, pz) + p;
pdz = dz;
dz = 2.0 * qmul(pz, pdz) + 1.0;
r = length(z);
dr = length(dz);
if (r > 8.0) break;
}
return 0.5 * log(r) * r / dr;
}
float de(vec3 p) {
return deQuaternionMandelbrot(vec4(p - vec3(0.5, 0.0, 0.0), 0.0));
}
参考
Quaternion Julia Set
http://glslsandbox.com/e#55839.0
ジュリア集合のクォータニオン版です。ジュリア集合は以下の漸化式で発散しない複素数(クォータニオン)$c$の集合になります。$a$は任意の値になります
\begin{align}
z_{0} &= c\\
z_{n+1} &= z_{n}^{2} + a
\end{align}
レイマーチングで必要なのは距離なので、マンデルブロ集合と同じようにDistance Estimationにより距離を求めます。
vec4 qmul(vec4 a, vec4 b) {
return vec4(
a.x * b.x - a.y * b.y - a.z * b.z - a.w * b.w,
a.x * b.y + a.y * b.x - a.z * b.w + a.w * b.z,
a.x * b.z + a.y * b.w + a.z * b.x - a.w * b.y,
a.x * b.w - a.y * b.z + a.z * b.y + a.w * b.x
);
}
#define ITERATIONS 16
float deQuaternionJuliaSet(vec4 p, vec4 c) {
vec4 z = p;
vec4 dz = vec4(1.0, 0.0, 0.0, 0.0);
vec4 pz, pdz;
float r = 0.0, dr = 1.0;
for (int i = 0; i < ITERATIONS; i++) {
pz = z;
z = qmul(pz, pz) + c;
pdz = dz;
dz = 2.0 * qmul(pz, pdz);
r = length(z);
dr = length(dz);
if (r > 4.0) break;
}
return 0.5 * log(r) * r / dr;
}
float de(vec3 p) {
return deQuaternionJuliaSet(vec4(p, 0.0), vec4(-1.0, 0.2, 0.0, 0.0));
}
参考
Mandelbulb
http://glslsandbox.com/e#55840.0
Mandelbulbです。マンデルブロ集合の漸化式の操作を極座標に変換してからおこなっているそうです。距離はマンデルブロ集合やジュリア集合と同じようにDistance Estimationにより求めています。詳しいことは理解していないので説明できません。
#define ITERATIONS 8
float deMandelbulb(vec3 p, float power) {
vec3 z = p;
float dr = 1.0;
float r;
for (int i = 0; i < ITERATIONS; i++) {
r = length(z);
if (r > 10.0) break;
float theta = acos(z.y / r);
float phi = atan(z.z, z.x);
dr = pow(r, power - 1.0) * power * dr + 1.0;
float zr = pow(r, power);
theta = theta * power;
phi = phi * power;
z = zr * vec3(sin(theta) * cos(phi), cos(theta), sin(theta) * sin(phi));
z += p;
}
return 0.5 * log(r) * r / dr;
}
float de(vec3 p) {
return deMandelbulb(p, 8.0);
}
参考
Mandelbox
http://glslsandbox.com/e#55841.0
box foldingとsphere foldingという折り畳みを行うフラクタルです。
#define foldingLimit 1.0
vec3 boxFold(vec3 z, float dz) {
return clamp(z, -foldingLimit, foldingLimit) * 2.0 - z;
}
void sphereFold(inout vec3 z, inout float dz, float minRadius, float fixedRadius) {
float m2 = minRadius * minRadius;
float f2 = fixedRadius * fixedRadius;
float r2 = dot(z, z);
if (r2 < m2) {
float temp = (f2 / m2);
z *= temp;
dz *= temp;
} else if (r2 < f2) {
float temp = (f2 / r2);
z *= temp;
dz *= temp;
}
}
// ref: http://blog.hvidtfeldts.net/index.php/2011/11/distance-estimated-3d-fractals-vi-the-mandelbox/
#define ITERATIONS 12
float deMandelbox(vec3 p, float scale, float minRadius, float fixedRadius) {
vec3 z = p;
float dr = 1.0;
for (int i = 0; i < ITERATIONS; i++) {
z = boxFold(z, dr);
sphereFold(z, dr, minRadius, fixedRadius);
z = scale * z + p;
dr = dr * abs(scale) + 1.0;
}
float r = length(z);
return r / abs(dr);
}
float de(vec3 p) {
return deMandelbox(p, 2.0, 0.5, 1.0);
}
参考
- Distance Estimated 3D Fractals (VI): The Mandelbox | Syntopia
- Mandelbox - Wikipedia
- 正解するカドの「カド」をレイマーチングでリアルタイム描画する | gam0022.net
Pseudo Kleinian
http://glslsandbox.com/e#55842.0
// ref: https://docs.google.com/presentation/d/1j4t4mcLw8F1PfqvKP3P8meMJ5dWfDXkfb9lc63qOFVM/edit#slide=id.g2460f5a976_0_0
#define ITERATIONS 8
float dePseudoKleinian(vec3 p) {
vec3 csize = vec3(0.90756, 0.92436, 0.90756);
float size = 1.0;
vec3 c = vec3(0.0);
float defactor = 1.0;
vec3 offset = vec3(0.0);
vec3 ap = p + 1.0;
for (int i = 0; i < ITERATIONS; i++) {
ap = p;
p = 2.0 * clamp(p, -csize, csize) - p;
float r2 = dot(p, p);
float k = max(size / r2, 1.0);
p *= k;
defactor *= k;
p += c;
}
float r = abs(0.5 * abs(p.y - offset.y) / defactor);
return r;
}
float de(vec3 p) {
return dePseudoKleinian(p);
}
参考
- [TDF2016] Indra's Bubbles
- CEDEC2017 TDF - Google スライド
- Distance Estimated 3D Fractals (Part VIII): Epilogue | Syntopia
- TokyoDemoFest2016に参加しました - 心鏡曼荼羅
終わりに
新しい距離関数があれば、随時追加していく予定です。もし面白い距離関数があれば教えていただけるとありがたいです。