Help us understand the problem. What is going on with this article?

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

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

Recursive Tetrahedron

recursive_tetrahedron.png
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

menger_sponge.png
メンガーのスポンジです。立方体から再帰的にクロス状の図形との論理差を取ることで穴を開けていきます。
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);
}

参考

Hexagonal Tiling

hexiagonal_tiling.png
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

quaternion_mandelbro_set.png
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

quaternion_julia_set.png
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

mandelbulb.png
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

mandelbox.png
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);
}

参考

Pseudo Kleinian

pesudo_kleinian.png
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);
}

参考

終わりに

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away