レイマーチングのための複雑な距離関数

レイマーチングで使える複雑めな距離関数を収集していきます。
球や直方体など単純な図形に関しては以下にまとまっているので、ここでは触れません。
[Inigo Quilez :: fractals, computer graphics, mathematics, shaders, demoscene and more](https://iquilezles.org/www/articles/distfunctions/distfunctions.htm)

## Recursive Tetrahedron

![recursive_tetrahedron.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/2c6ebac7-f941-d51e-2926-64a187f6cd5a.png)
http://glslsandbox.com/e#55834.0
四面体のIFSによる繰り返しです。最終的に四面体の頂点に対して距離を計算しています。

```glsl
#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);
}
```

### 参考

* [Distance Estimated 3D Fractals \(III\): Folding Space \| Syntopia](http://blog.hvidtfeldts.net/index.php/2011/08/distance-estimated-3d-fractals-iii-folding-space/)

## Menger Sponge

![menger_sponge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/163a8fbc-b915-14a4-3490-07c2fe1cfc07.png)
メンガーのスポンジです。立方体から再帰的にクロス状の図形との論理差を取ることで穴を開けていきます。
http://glslsandbox.com/e#55835.0

```glsl
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

```glsl
// 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](https://qiita.com/aa_debdeb/items/03609566ea322a4b764f)
* [Inigo Quilez :: fractals, computer graphics, mathematics, shaders, demoscene and more](https://www.iquilezles.org/www/articles/menger/menger.htm)
* [Unity×レイマーチングによる映像制作の実践手法 \| gam0022\.net](https://gam0022.net/blog/2019/06/25/unity-raymarching/)

## Hexagonal Tiling

![hexiagonal_tiling.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/463a5cf1-9e93-0f92-02e8-b3efc4ec8c8a.png)
http://glslsandbox.com/e#55837.0
六角形のタイリングです。六角形ではグリッド状にきれいに分割することができないので、単純に`mod`を用いた繰り返しを使用してもうまくいきません。そのため、六角形のタイリングを2つに分けてそれぞれ距離を計算し、論理和で結合しています。

```glsl
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

![quaternion_mandelbro_set.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/254de56a-96b1-6f5f-2538-6901e45e998a.png)
http://glslsandbox.com/e#55838.0
マンデルブロ集合のクォータニオン版です。マンデルブロ集合は以下の漸化式で発散しない複素数(クォータニオン)$c$の集合になります。

```math
\begin{align}
z_{0} &= 0\\
z_{n+1} &= z_{n}^{2} + c
\end{align}
```

この漸化式では$c$で発散するかどうかしかわかりません。レイマーチングでは物体表面までの距離が必要になりますが、Distance Estimationという手法により距離の下限値は以下の式で求まります。

```math
DE = 0.5r(\ln{r}) / dr
```

ここで$r$は`escape time length`、$dr$は`length of running derivative`になり、以下により求まります。
(私も詳しく説明できないので、詳細に知りたい方は[この記事](http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-iv-the-holy-grail/)を読んでください)

```math
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)|
```

これを愚直に実装したものが以下のようになります。漸化式の積がクォータニオンのものになっていることに注意してください。

```glsl

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));
}
```


### 参考

* [Distance Estimated 3D Fractals \(IV\): The Holy Grail \| Syntopia](http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-iv-the-holy-grail/)

## Quaternion Julia Set

![quaternion_julia_set.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/4203d21b-349e-9291-9d7f-726391ff42ed.png)
http://glslsandbox.com/e#55839.0
ジュリア集合のクォータニオン版です。ジュリア集合は以下の漸化式で発散しない複素数(クォータニオン)$c$の集合になります。$a$は任意の値になります

```math
\begin{align}
z_{0} &= c\\
z_{n+1} &= z_{n}^{2} + a
\end{align}
```

レイマーチングで必要なのは距離なので、マンデルブロ集合と同じようにDistance Estimationにより距離を求めます。

```glsl
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));
}
```

### 参考

* [Distance Estimated 3D Fractals \(IV\): The Holy Grail \| Syntopia](http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-iv-the-holy-grail/)
* [Quaternion Julia Fractals](http://paulbourke.net/fractals/quatjulia/)



## Mandelbulb

![mandelbulb.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/88d4f38b-0c17-4010-53db-f2901ea04c0e.png)
http://glslsandbox.com/e#55840.0
Mandelbulbです。マンデルブロ集合の漸化式の操作を極座標に変換してからおこなっているそうです。距離はマンデルブロ集合やジュリア集合と同じようにDistance Estimationにより求めています。詳しいことは理解していないので説明できません。

```glsl
#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);
}
```

### 参考

* [Distance Estimated 3D Fractals \(V\): The Mandelbulb & Different DE Approximations \| Syntopia](http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations/)

## Mandelbox

![mandelbox.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/0595a6f5-4fa8-1bbf-435c-11ebfbff8ec8.png)
http://glslsandbox.com/e#55841.0
box foldingとsphere foldingという折り畳みを行うフラクタルです。

```glsl
#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](http://blog.hvidtfeldts.net/index.php/2011/11/distance-estimated-3d-fractals-vi-the-mandelbox/)
* [Mandelbox \- Wikipedia](https://en.wikipedia.org/wiki/Mandelbox)
* [正解するカドの「カド」をレイマーチングでリアルタイム描画する \| gam0022\.net](https://gam0022.net/blog/2017/06/30/raymarching-kado/)

## Pseudo Kleinian

![pesudo_kleinian.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/b5c264f0-6a23-0ad4-8f8f-4889d938588b.png)
http://glslsandbox.com/e#55842.0

```glsl
// 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](https://www.shadertoy.com/view/XsGGWG)
* [CEDEC2017 TDF \- Google スライド](https://docs.google.com/presentation/d/1j4t4mcLw8F1PfqvKP3P8meMJ5dWfDXkfb9lc63qOFVM/edit#slide=id.g2460f5a976_0_0)
* [Distance Estimated 3D Fractals \(Part VIII\): Epilogue \| Syntopia](http://blog.hvidtfeldts.net/index.php/2012/05/distance-estimated-3d-fractals-part-viii-epilogue/)
* [TokyoDemoFest2016に参加しました \- 心鏡曼荼羅](http://soma.hatenablog.jp/entry/2016/02/22/220100)


## 終わりに

新しい距離関数があれば、随時追加していく予定です。もし面白い距離関数があれば教えていただけるとありがたいです。