GLSLレイマーチング研究_距離関数について勉強してみた05(Torusの関数をいじる(距離関数導出編))

参考記事:[modeling with distance functions](http://iquilezles.org/www/articles/distfunctions/distfunctions.htm)

三記事にも及ぶ、Box関数が終わりようやく、次のTorusに行けます。

![Torus_base.PNG](https://qiita-image-store.s3.amazonaws.com/0/102634/1c348116-27f9-279c-5b5c-2ddeb944fcc3.png)

## 今回ですが、
距離関数の導出に分量がかかってしまったので、記事2回にわけます。

## Torusってなによ。
簡単にいうとドーナツの形です。

数学やっていると、幾何でよく見ます…

よく、異世界の形で例えられます。

例えば、ドラクエの世界は、Torusとかいわれますね。

説明面倒なので、この記事見てね.

[ドラゴンクエストのトーラス世界](http://d.hatena.ne.jp/LM-7/20071219/1199055762)

## いきなり、距離関数はきつい…
今回は、少し式が複雑なので、苦手な人もいると思いますが…

数式から書いています…（数学アレルギーの人は、読み飛ばしても…）

## Torusの数式

```math
S(u,v) = 
\begin{pmatrix}
(R + r \cos u) \cos v \\
(R + r \cos u) \sin v \\
r \sin u
\end{pmatrix}
, (0 \leq u, v \leq 2 \pi), (0 < r < R)
```

### 変数の補足

```math
\begin{align}
u&:円環の断面内の角度 \\
v&:円環の中心方向の角度 \\
r&:小半径 \\
R&:大半径 \\
\end{align}
```

てな感じです.

### Torusの距離関数の導出
今回,久々にガチで計算してTorusの距離関数を導出してみました…（興味がある人は見てね。）

```math
x = (R + r \cos u) \cos v, y = (R + r \cos u) \sin v, z = r \sin u
```

```math
\begin{align}
x^2+y^2 &= (R + r \cos u)^2 \cos^2 u + (R + r \cos u)^2 \sin^2 v \\
&= (R + r \cos u)^2( \cos^2 u + \sin^2 v) \\
&= (R + r \cos u)^2 \\
&= R^2 + 2Rr \cos u + r^2 \cos^2 u
\end{align}
```

ここで,

```math
r \cos u = \sqrt{r^2+z^2}
```
より

```math
\begin{align}
x^2+y^2 &= R^2 + 2R (r \cos u) + (r \cos u)^2 \\
&= R^2 + 2R \sqrt{r^2+z^2} + (r^2+z^2)
\end{align}
```

まとめると

```math
\begin{align}
(x^2+y^2+z^2) - R^2 - r^2 &= 2R \sqrt{r^2+z^2} \\
\end{align}
```

両辺を二乗

```math
\begin{align}
((x^2+y^2+z^2-R^2) - r^2 )^2 &= 4R^2 (r^2+z^2)
\end{align}
```

`r`についてまとめると（ここが割とダルい…）

```math
\begin{align}
r^4 -(x^2+y^2+z^2+R^2) r^2 + (x^2+y^2+z^2+R^2)^2 - 2 R^2 (x^2+y^2+z^2-R^2) + R^4 = 0
\end{align}
```

`r^2` について[二次方程式の解の公式](http://w3e.kanazawa-it.ac.jp/math/category/kansuu/2jihouteisiki/henkan-tex.cgi?target=/math/category/kansuu/2jihouteisiki/kainokousiki.html)で解くと

```math
\begin{align}
r^2 &= (x^2+y^2+z^2+R^2) \pm \sqrt{(x^2+y^2+z^2+R^2)^2 - (x^2+y^2+z^2)^2 + 2 R^2 (x^2+y^2+z^2) - R^4} \\
&= (x^2+y^2+z^2+R^2) \pm \sqrt{4 R^2 (x^2+y^2)} \\
&= (x^2+y^2+z^2+R^2) \pm 2R \sqrt{x^2+y^2} \\
ここで, (0 < r < R) より \\
r^2 &= (x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}
\end{align}
```

最後に,

```math
((x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}) -r^2 = 0 \\
(\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} -r)(\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} +r) = 0 \\
ここで、(\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} +r)は、rが、視界に入らないので,不適
```
```math
\begin{align}
(\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} -r) = 0 
\end{align}
```

となり、距離関数は、

```math
\begin{align}
\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} -r
\end{align}
```
となります.

ちなみに、`length` を使った関数の導出も書いておきます.

```math
\begin{align}
\sqrt{(x^2+y^2+z^2+R^2) - 2R \sqrt{x^2+y^2}} -r &= \sqrt{(x^2+y^2) - 2R \sqrt{x^2+y^2} + R^2 +z^2} -r \\
&= \sqrt{((\sqrt{x^2+y^2})^2 - 2R \sqrt{x^2+y^2} + R^2) + z^2} -r \\
&= \sqrt{(\sqrt{x^2+y^2}) - R)^2 + z^2} -r \\
\end{align}
```

ルートを `length` に書き直すと得られます.

### 書き方の説明

Torusは、xz平面に半径 `r` の円があり、その円をy軸で半径`R`で回転させた図形です…

## さっそく、コードにすると

```c++
float r = 【小半径】;
float R = 【大半径】;
return sqrt(p.x*p.x+p.y*p.y+p.z*p.z+R*R - 2.0 * R * sqrt(p.x*p.x+p.y*p.y) ) - r;
```

`length`を使って書き直すと

```c++
vec2 t = vec2(【大半径】, 【小半径】);
vec2 q = vec2(length(p.xz)-t.x,p.y);
return length(q)-t.y;
```

すっきり書けますね.

![Torus_base02.PNG](https://qiita-image-store.s3.amazonaws.com/0/102634/0b6aa22c-7a15-6a2e-325f-4ce668c1812b.png)

コード

```c++
// ============================================================================
// Torus function
// ============================================================================

precision mediump float;
uniform vec2  resolution;    // resolution (512.0, 512.0)
uniform vec2  mouse;         // mouse      (-1.0 ~ 1.0)
uniform float time;          // time       (1second == 1.0)
uniform sampler2D prevScene; // previous scene texture

// Torusの距離関数
float dTorus(vec3 p){
    // vec2 t = vec2(1.0, 0.3);
    // vec2 q = vec2(length(p.xz)-t.x,p.y);
    // return length(q)-t.y;
    
    float r = 0.3;
    float R = 1.0;
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z+R*R - 2.0 * R * sqrt(p.x*p.x+p.y*p.y) ) - r;

}

// 距離関数を呼び出すハブ関数
float distanceHub(vec3 p){
    return dTorus(p);
    // return max(dCube(p), length(p)-1.0*abs(sin(time)));
}

// 法線を生成する
vec3 genNormal(vec3 p){
    float d = 0.001;
    return normalize(vec3(
        distanceHub(p + vec3(  d, 0.0, 0.0)) - distanceHub(p + vec3( -d, 0.0, 0.0)),
        distanceHub(p + vec3(0.0,   d, 0.0)) - distanceHub(p + vec3(0.0,  -d, 0.0)),
        distanceHub(p + vec3(0.0, 0.0,   d)) - distanceHub(p + vec3(0.0, 0.0,  -d))
    ));
}

void main(){
    // スクリーンスペースを考慮して座標を正規化する
    vec2 p = (gl_FragCoord.xy * 2.0 - resolution) / min(resolution.x, resolution.y);

    // カメラを定義する
    vec3 cPos = vec3(0.0,  1.0,  2.5); // カメラの位置
    vec3 cDir = vec3(0.0,  -0.5, -1.0); // カメラの向き(視線)
    vec3 cUp  = vec3(0.0,  1.0,  0.0); // カメラの上方向
    vec3 cSide = cross(cDir, cUp);     // 外積を使って横方向を算出
    float targetDepth = 1.0;           // フォーカスする深度

    // カメラの情報からレイを定義する
    vec3 ray = normalize(cSide * p.x + cUp * p.y + cDir * targetDepth);

    // マーチングループを組む
    float dist = 0.0;  // レイとオブジェクト間の最短距離
    float rLen = 0.0;  // レイに継ぎ足す長さ
    vec3  rPos = cPos; // レイの先端位置(初期位置)
    for(int i = 0; i < 32; ++i){
        dist = distanceHub(rPos);
        rLen += dist;
        rPos = cPos + ray * rLen;
    }

    // レイとオブジェクトの距離を確認
    if(abs(dist) < 0.001){
        // 法線を算出
        vec3 normal = genNormal(rPos);

        // ライトベクトルの定義（マウスの影響を受けるように）
        // vec3 light = normalize(vec3(mouse + 1.0, 1.0));
        vec3 light = normalize(vec3(1.0, 1.0, 1.0));

        // ライトベクトルとの内積を取る
        float diff = max(dot(normal, light), 0.1);

        // diffuse を出力する
        gl_FragColor = vec4(vec3(diff), 1.0);
    }else{
        // 衝突しなかった場合はそのまま黒
        gl_FragColor = vec4(vec3(0.0, 0.0, 0.0), 1.0);
    }
}
```