[コードリーディング vol.2] レイマーチングによる雲表現を読み解く

## 概要

以前に、「[[コードリーディング] レイマーチングによる波表現を読み解く](https://qiita.com/edo_m18/items/a575606a60b21f0d2c57)」というタイトルで波表現をするレイマーチングのコードを読み解きました。
（そのときの動画はこれ↓）
![レイマーチで描く海](https://qiita-image-store.s3.amazonaws.com/0/5616/e44cda7e-83c0-69dc-8373-1bf7bcdfffec.gif)

今回は[こちらのShadertoyの作品](https://www.shadertoy.com/view/lss3zr)をコードリーディングしたいと思います。
ちなみにこんな感じの絵が出ます↓
（下記は、上の作品にコメントつけたり分かりやすくしたりしたもののキャプチャです）

![demo.gif](https://qiita-image-store.s3.amazonaws.com/0/5616/e17538c5-4f2b-2982-7374-501c77726bfb.gif)



こちらの作品を、自分が分かる範囲でコメントを付けつつ、変数などもできるだけ意味付けしながら解説・読解していこうと思います。

なお、[自分がコメントを付けたバージョン](https://www.shadertoy.com/view/WslGWl)はこちらにポストしてあります。

また読解に際して、こちらのブログ記事（[Creating a Volumetric Ray Marcher](https://shaderbits.com/blog/creating-volumetric-ray-marcher)）も参考にさせていただきました。

いちおう翻訳とかも全文やってみたのですが、書いた方が許可してくれたら掲載しようと思います。


許可されたので掲載しました。

- [Cloud Raymachingの翻訳記事](https://qiita.com/edo_m18/items/876f2857e67e26a053d6)



## コード全文

まずはコード全文を掲載します。
見てもらうと分かりますが、コード量はそれほど多くありません。

```glsl
// ------------------------------------------
//
// This post cloned from "https://www.shadertoy.com/view/lss3zr"
//
// I also refer this blog post below.
// https://shaderbits.com/blog/creating-volumetric-ray-marcher
//
// This post is to learn how to cloud raymarching is working.
//
// ------------------------------------------

#define USE_LIGHT 0

mat3 m = mat3( 0.00,  0.80,  0.60,
              -0.80,  0.36, -0.48,
              -0.60, -0.48,  0.64);

float hash(float n)
{
    return fract(sin(n) * 43758.5453);
}

///
/// Noise function
///
float noise(in vec3 x)
{
    vec3 p = floor(x);
    vec3 f = fract(x);
    
    f = f * f * (3.0 - 2.0 * f);
    
    float n = p.x + p.y * 57.0 + 113.0 * p.z;
    
    float res = mix(mix(mix(hash(n +   0.0), hash(n +   1.0), f.x),
                        mix(hash(n +  57.0), hash(n +  58.0), f.x), f.y),
                    mix(mix(hash(n + 113.0), hash(n + 114.0), f.x),
                        mix(hash(n + 170.0), hash(n + 171.0), f.x), f.y), f.z);
    return res;
}

///
/// Fractal Brownian motion.
///
/// Refer to:
/// EN: https://thebookofshaders.com/13/
/// JP: https://thebookofshaders.com/13/?lan=jp
///
float fbm(vec3 p)
{
    float f;
    f  = 0.5000 * noise(p); p = m * p * 2.02;
    f += 0.2500 * noise(p); p = m * p * 2.03;
    f += 0.1250 * noise(p);
    return f;
}

//////////////////////////////////////////////////

///
/// Like sphere distance function.
///
/// But this function return inverse value.
/// Normal dist function is like below.
/// 
/// return length(pos) - 0.1;
///
/// Because this function is used for density.
///
float scene(in vec3 pos)
{
    return 0.1 - length(pos) * 0.05 + fbm(pos * 0.3);
}

///
/// Get normal of the cloud.
///
vec3 getNormal(in vec3 p)
{
    const float e = 0.01;
    return normalize(vec3(scene(vec3(p.x + e, p.y, p.z)) - scene(vec3(p.x - e, p.y, p.z)),
                          scene(vec3(p.x, p.y + e, p.z)) - scene(vec3(p.x, p.y - e, p.z)),
                          scene(vec3(p.x, p.y, p.z + e)) - scene(vec3(p.x, p.y, p.z - e))));
}

///
/// Create a camera pose control matrix.
///
mat3 camera(vec3 ro, vec3 ta)
{
    vec3 cw = normalize(ta - ro);
    vec3 cp = vec3(0.0, 1.0, 0.0);
    vec3 cu = cross(cw, cp);
    vec3 cv = cross(cu, cw);
    return mat3(cu, cv, cw);
}

///
/// Main function.
///
void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = (fragCoord.xy * 2.0 - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    vec2 mo = vec2(iTime * 0.1, cos(iTime * 0.25) * 3.0);
    
    // Camera
    float camDist = 25.0;
    
    // target
    vec3 ta = vec3(0.0, 1.0, 0.0);
    
    // Ray origin
    //vec3 ori = vec3(sin(iTime) * camDist, 0, cos(iTime) * camDist);
    vec3 ro = camDist * normalize(vec3(cos(2.75 - 3.0 * mo.x), 0.7 - 1.0 * (mo.y - 1.0), sin(2.75 - 3.0 * mo.x)));
    
    float targetDepth = 1.3;
    
    // Camera pose.
    mat3 c = camera(ro, ta);
    vec3 dir = c * normalize(vec3(uv, targetDepth));
    
    // For raymarching const values.
    const int sampleCount = 64;
    const int sampleLightCount = 6;
    const float eps = 0.01;
    
    // Raymarching step settings.
    float zMax = 40.0;
    float zstep = zMax / float(sampleCount);
    
    float zMaxl = 20.0;
    float zstepl = zMaxl / float(sampleLightCount);
    
    // Easy access to the ray origin
    vec3 p = ro;
    
    // Transmittance
    float T = 1.0;
    
    // Substantially transparency parameter.
    float absorption = 100.0;
    
    // Light Direction
    vec3 sun_direction = normalize(vec3(1.0, 0.0, 0.0));
    
    // Result of culcration
    vec4 color = vec4(0.0);
    
    for (int i = 0; i < sampleCount; i++)
    {
        // Using distance function for density.
        // So the function not normal value.
        // Please check it out on the function comment.
        float density = scene(p);
        
        // The density over 0.0 then start cloud ray marching.
        // Why? because the function will return negative value normally.
        // But if ray is into the cloud, the function will return positive value.
        if (density > 0.0)
        {
            // Let's start cloud ray marching!
            
            // why density sub by sampleCount?
            // This mean integral for each sampling points.
            float tmp = density / float(sampleCount);
            
            T *= 1.0 - (tmp * absorption);
            
            // Return if transmittance under 0.01. 
            // Because the ray is almost absorbed.
            if (T <= 0.01)
            {
                break;
            }
            
            #if USE_LIGHT == 1
            // Light scattering
            
            // Transmittance for Light
            float Tl = 1.0;
            
            // Start light scattering with raymarching.
            
            // Raymarching position for the light.
            vec3 lp = p;
            
            // Iteration of sampling light.
            for (int j = 0; j < sampleLightCount; j++)
            {
                float densityLight = scene(lp);
                
                // If densityLight is over 0.0, the ray is stil in the cloud.
                if (densityLight > 0.0)
                {
                    float tmpl = densityLight / float(sampleCount);
                    Tl *= 1.0 - (tmpl * absorption);
                }
                
                if (Tl <= 0.01)
                {
                    break;
                }
                
                // Step to next position.
                lp += sun_direction * zstepl;
            }
            #endif
            
            // Add ambient + light scattering color
            float opaity = 50.0;
            float k = opaity * tmp * T;
            vec4 cloudColor = vec4(1.0);
            vec4 col1 = cloudColor * k;
            
            #if USE_LIGHT == 1
            float opacityl = 30.0;
            float kl = opacityl * tmp * T * Tl;
            vec4 lightColor = vec4(1.0, 0.7, 0.9, 1.0);
            vec4 col2 = lightColor * kl;
            #else
            vec4 col2 = vec4(0.0);
            #endif
            
            color += col1 + col2;
        }
        
        p += dir * zstep;
    }
    
    vec3 bg = mix(vec3(0.3, 0.1, 0.8), vec3(0.7, 0.7, 1.0), 1.0 - (uv.y + 1.0) * 0.5);
    color.rgb += bg;
    
    fragColor = color;
}
```


さて、ではこれをベースにひとつずつ見ていきましょう。


## ノイズ関数とFBM

冒頭で定義されている関数を見てみましょう。

```glsl
mat3 m = mat3( 0.00,  0.80,  0.60,
              -0.80,  0.36, -0.48,
              -0.60, -0.48,  0.64);

float hash(float n)
{
    return fract(sin(n) * 43758.5453);
}

///
/// Noise function
///
float noise(in vec3 x)
{
    vec3 p = floor(x);
    vec3 f = fract(x);
    
    f = f * f * (3.0 - 2.0 * f);
    
    float n = p.x + p.y * 57.0 + 113.0 * p.z;
    
    float res = mix(mix(mix(hash(n +   0.0), hash(n +   1.0), f.x),
                        mix(hash(n +  57.0), hash(n +  58.0), f.x), f.y),
                    mix(mix(hash(n + 113.0), hash(n + 114.0), f.x),
                        mix(hash(n + 170.0), hash(n + 171.0), f.x), f.y), f.z);
    return res;
}

///
/// Fractal Brownian motion.
///
/// Refer to:
/// EN: https://thebookofshaders.com/13/
/// JP: https://thebookofshaders.com/13/?lan=jp
///
float fbm(vec3 p)
{
    float f;
    f  = 0.5000 * noise(p); p = m * p * 2.02;
    f += 0.2500 * noise(p); p = m * p * 2.03;
    f += 0.1250 * noise(p);
    return f;
}
```

### FBMで使われるオフセット行列

冒頭で定義されている行列（`m`）はおそらくオフセット行列でしょう。
ノイズに引数として与えられる点を少しずつ回転させる意味がありそうです。

実際に2x2行列にして点の移動を見てみると、適当にオフセットさせるような動きをします。
desmosで実際に点を動かしてみると、行列を適用した別の点がどう動くかを確認することができます。

[![](https://qiita-image-store.s3.amazonaws.com/0/5616/1455f827-eb3c-139b-6a30-d3dbc1c0b606.png)](https://www.desmos.com/calculator/e0tlzb36sy)[https://www.desmos.com/calculator/e0tlzb36sy](https://www.desmos.com/calculator/e0tlzb36sy)


### hash関数とnoise関数

次に記述されているのは`hash`関数と`noise`関数です。
これらに関しては前回のコードリーディングの記事で似た部分について言及しているのでそちらを見てもらうといいでしょう。

- [[コードリーディング] レイマーチングによる波表現を読み解く](https://qiita.com/edo_m18/items/a575606a60b21f0d2c57#hash%E3%81%A8%E5%90%8D%E3%81%AE%E3%81%A4%E3%81%8F%E3%83%AF%E3%82%B1)


ざっくりとだけ説明しておくと、`hash`関数は引数が同じなら同じ値を返す関数です。
ここで使いたいものは「疑似乱数」でしょう。そして疑似乱数は`hash`関数同様、引数と種が同じなら同じ値を返す関数です。
つまりは「疑似乱数」を得ている、と見ていいでしょう。


### FBM

そして冒頭での最後の関数は`fbm`関数です。
これは「**Fractal Brownian motion**」の略で、非整数ブラウン運動と訳されます。
これについては以前記事に書いたのでそちらを参照ください。

- [フラクタルブラウン運動とドメインワープ](https://qiita.com/edo_m18/items/e4d7a084cdbbfdc7863c)

要は、雲のような波のうねり表現を加えてくれるもの、ですね。


### 距離関数の亜種

さて、次に見るのは`scene`関数です。
これは距離関数の亜種として利用されています。
`scene`関数は以下のように定義されています。

```glsl:距離関数？
float scene(in vec3 pos)
{
    return 0.1 - length(pos) * 0.05 + fbm(pos * 0.3);
}
```

ちょっと見慣れない形ですね。
ただ、ここで行っている処理をざっくりいうと、通常の球体の距離関数`fbm`関数でのノイズを混ぜたもの、という感じです。


#### ただし密度を計算している

少し注意してもらいたいのは、距離関数は本来、「とある点から該当オブジェクトまでの距離」を返す関数であるということです。
そして関数をよく見てみると、主な距離関数部分を「減算」する形になっていることに気づきます。

```glsl:一部抜粋
return 0.1 - length(pos);
```

さて、普通の球体の距離関数は以下です。

```glsl:球体の距離関数
return length(pos) - 1.0;
```

計算が逆転していることが分かりますね。
そして、離れた位置からの計算をした場合（つまり上記関数で言えば`0.1`よりも離れた場所の場合）、計算結果はマイナスとなってしまいます。

これは、この関数が通常の距離関数ではく「密度」を返しているからです。
つまり、レイが球体の中に入った（`0.1`以下になった）ら、その値はプラスに転じます。
そしてレイが進むに連れてその値が変化し、仮に中心に向かっていくレイの場合はその値はマーチを繰り返すたびに大きくなり、さらに中心を抜けるとまた徐々に小さくなっていきます。

図にすると以下のような感じですね。

![密度の遷移](https://qiita-image-store.s3.amazonaws.com/0/5616/e9ddd0a7-6090-8285-53e8-f48e5a40552a.jpeg)

**step size**分ループするたびに進み、近づいていくにつれて関数が返す値は徐々に`0`に近づき、オブジェクトの中に入ったらプラスに転じる様子を示しています。


#### レイは「一定距離ずつ」進む

さて、もうひとつ注意点があります。
それは、通常のレイマーチングでは「一番近いオブジェクト」への距離分、レイが移動していきます。
しかし今回のレイマーチングではそうではなく、あらかじめ規定した距離（＝一定距離）ずつ進んでいくことになります。

これは、不透明オブジェクトを描き出すのが目的ではなく、あくまで半透明である「雲」をレンダリングするため、一定の距離ごとにサンプリング点を取り、その位置からライト方向への密度をサンプリングしそれを累積することで雲の絵を描き出す、という手法のためです。

そのため、オブジェクトごとの最短距離、という判定ではなく一定距離ずつレイを進める必要がある、というわけなんですね。


### main関数

さてではいよいよ、main関数の中身を見ていきましょう。
といっても、冒頭は通常のレイマーチと同様に、テクセルの位置を元にカメラ姿勢を適用してレイの方向を決める処理になります。
このあたりは以前書いたレイマーチング入門などにも書いているので詳細はそちらをご覧ください。

- [[GLSL] レイマーチング入門 vol.1](https://qiita.com/edo_m18/items/034665d42c562da88cb6)

今回は雲のレンダリング処理部分に焦点を当てて書いていきたいと思います。


```glsl:変数セットアップ
// For raymarching const values.
const int sampleCount = 64;
const int sampleLightCount = 6;
const float eps = 0.01;

// Raymarching step settings.
float zMax = 40.0;
float zstep = zMax / float(sampleCount);

float zMaxl = 20.0;
float zstepl = zMaxl / float(sampleLightCount);

// レイの視点（Ray Origin）を簡易アクセスするために`p`を定義
vec3 p = ro;

// Transmittance
float T = 1.0;

// 実質的には透明度
float absorption = 100.0;

// ライト方向
vec3 sun_direction = normalize(vec3(1.0, 0.0, 0.0));

// 計算結果用
vec4 color = vec4(0.0);
```

まず最初は変数のセットアップです。


#### サンプル回数を定義

`sampleCount`はレイマーチング自体のサンプルする点の数、`sampleLightCount`は、左記のサンプル点からライトへ向かうベクトル中のサンプルの数となっています。

ちなみにこのサンプル点のイメージは以下になります。

![サンプルリングイメージ](https://qiita-image-store.s3.amazonaws.com/0/5616/01ba0960-6d71-891e-6f40-3b0bf8a762d3.jpeg)


まず、レイが前述の通り一定距離（**step size**）ずつ進んでいきます。
そしてレイがオブジェクト内に入っていたら、ライト方向に対してもサンプリングを行います。
ちなみにこれは最大、N * Mという計算量になり、とても高価な処理になります。
（N回のサンプル点に対して、ベクトル方向にもM回サンプルするため、最大でN * M回のサンプルが必要となります）

今回の例では64回のレイマーチングによるサンプルと、ライト方向への6回のサンプルを実行しているため、最大で384回のサンプリングが発生する計算になります。
距離関数が単純であればそこまで高負荷はなりませんが、複雑になればなるほどその影響が大きくなっていきます。


#### サンプル距離を計算

サンプル回数と最大サンプリング距離が決まれば自ずと一回のマーチで進むレイの距離が計算できます。
それを行っているのが以下の計算部分です。

```glsl:レイの進む距離を事前計算
// Raymarching step settings.
float zMax = 40.0;
float zstep = zMax / float(sampleCount);

float zMaxl = 20.0;
float zstepl = zMaxl / float(sampleLightCount);
```

上が通常のレイマーチング、下がライトへ向かうサンプルの距離です。
単純に、定めた距離をサンプル回数で割っているだけですね。

#### 透過率・吸収率

さて、最後は以下の変数です。
（まだ説明していない変数もありますが、名称から分かる or コメントで説明できる範囲のものはコメントにしてあるのでそちらを見てください）

```glsl:透過関連の設定
// Transmittance
float T = 1.0;

// Substantially transparency parameter.
float absorption = 100.0;
```

ここで`T`は**Transmittance**、つまり「透過率」ですね。透過率を計算結果に掛けることで徐々に不透明な感じを出していくための変数となります。

次の`absorption`は、英単語だけで言えば**吸収**という意味なので、光が吸収されたことを計算するために用いられるものだと思います。
つまり「**どれだけ光を吸収するか**」の値ですね。

計算部分を見てみると以下のようになっています。

```glsl:absorptionの計算
T *= 1.0 - (tmp * absorption);
```

となっており、「密度 * absorption」を「1.0」から減算したものを透過率`T`に掛けています。

吸収率`absorption`の値が大きければ大きいほど光が吸収されるため、レイマーチング的には早々にリターンすることになります。
なので値が大きいほど雲は「透明に近く」なります。
実質的には「透明度（**Transparency**）」ということになりますね。


### レイマーチング

さて、ではいよいよレイマーチング部分を見ていきましょう。

通常のレイマーチングと同様、指定した回数`for`文によってループさせます。
冒頭では距離関数を利用して「密度」を求めます。

```glsl:距離関数の亜種を用いて密度を計算
float density = scene(p);
```

そしてそれがプラスに転じたらサンプリング開始です。

```glsl:密度の判定
if (density > 0.0) { ... }
```

距離関数の亜種のところでも書きましたが、この関数は「密度」を返します。
つまり「雲の中」にレイがある場合は密度がある状態、「雲の外」にレイがある場合は密度がマイナスとすることで、うまく「雲の中」でだけサンプリングできるようにしている、というわけなんですね。
なので、プラスに転じた段階でサンプリングを開始しているわけです。

ちなみに、密度がマイナスの間は一定の距離ずつレイを進めていきます。

```glsl:レイを進める
p += dir * zstep;
```

`dir`がレイの方向、`zstep`が事前に計算した1stepで進む距離です。
それを加算していくことでレイを進めています。


#### 透過率計算

```glsl:透過率の計算
float tmp = density / float(sampleCount);

T *= 1.0 - (tmp * absorption);
```

透過率`T`は、各サンプル単位に正規化（`sampleCount`で除算）した値を使って計算します。
また上で書いたように、`absorption`を利用してどれくらい透過していくかも計算します。

そして透過率がしきい値（今回の例では`0.01`）を下回った場合はレイマーチングを終了します。

```glsl
// Return if transmittance under 0.01. 
// Because the ray is almost absorbed.
if (T <= 0.01)
{
    break;
}
```

### 色の計算

実は密度だけの計算でいいならこれで計算は終わりです。
なので、`#define`を使ってライトの有無を切り替えられるようにしてあります。

試しに`USE_LIGHT`を`1`にしてみると色が加算され、より雲っぽくなります。

まずは光の計算部分をいったんスキップして、色を計算する部分を見てみましょう。


```glsl:色の計算
float opaity = 50.0;
float k = opaity * tmp * T;
vec4 cloudColor = vec4(1.0);
vec4 col1 = cloudColor * k;

// ... 中略 ...

color += col1;
```

`cloudColor`が雲自体の色です。上の例では「白`(1.0, 1.0, 1.0)`」ですね。

そして係数の計算は（実質的な）密度`tmp`に、透過率`T`を掛け、さらにそれを`50.0`倍して、それを色に乗算することで求めています。
この`50.0`は実質的に不透明度（opaity）ということになりますね。

なのでこの値を小さくすると雲の色が薄くなります。

ここまでで、冒頭の雲の絵が描かれることになります。


### 光の計算を加える

さて、では光を加える計算を見てみましょう。以下の部分ですね。


```glsl:光の計算
// Light scattering

// Transmittance for Light
float Tl = 1.0;

// Start light scattering with raymarching.

// Raymarching position for the light.
vec3 lp = p;

// Iteration of sampling light.
for (int j = 0; j < sampleLightCount; j++)
{
    float densityLight = scene(lp);

    // If densityLight is over 0.0, the ray is stil in the cloud.
    if (densityLight > 0.0)
    {
        float tmpl = densityLight / float(sampleCount);
        Tl *= 1.0 - (tmpl * absorption);
    }

    if (Tl <= 0.01)
    {
        break;
    }

    // Step to next position.
    lp += sun_direction * zstepl;
}
```

冒頭の変数たちはライト計算に使う透過率とサンプリングに使うレイの開始位置です。
（レイマーチング本体の位置を上書かないようコピーしています）

そして通常のレイマーチと同様に、光方向に向かってサンプリングを行っていきます。
レイマーチの計算はほぼ、本体のレイマーチングと同様ですね。方向がライト方向になっているのが異なるくらいでしょう。

そして本体と同様、透過率がしきい値以下になったらループを抜けています。

### 光の色を加算する

最後に、上記までで計算した`Tl`を利用して光の色を加算します。

```glsl:光の色の計算
float opacityl = 30.0;
float kl = opacityl * tmp * T * Tl;
vec4 lightColor = vec4(1.0, 0.7, 0.9, 1.0);
vec4 col2 = lightColor * kl;
```

計算式自体は密度の計算で利用したものとほぼ同じですね。
違いは、係数の計算に`T`も利用して、`T * Tl`となっている点くらいです。


### 色の計算

以上でコード全容を見てきました。
最後の色の計算部分について考えてみます。

計算式は、言葉で表すと

「**不透明度 * 密度 * 透過率 * 色**」

という形になっています。

透過率自体は以下の式になっていました。

```glsl:透過率の計算
float tmp = density / float(sampleCount);
T *= 1.0 - (tmp * absorption);
```

1.0から密度に吸収率を掛けたものを引いたものを累積（`*=`）したものですね。
おそらく、0.0〜1.0の間であることを期待している式だと思います。
（実際、試しに`clamp(tmp * absorption, 0.0, 1.0)`を追加してみましたが結果は変わりませんでした）

つまり、透過率は一方的に減っていく（少数値のみが掛けられる）ことになります。
そしてその累積した透過率がしきい値（`0.01`）を下回ったら、光がそれ以上進まないことを意味するのでループを終了しているわけですね。


## まとめ

以上で雲のレイマーチングによるレンダリング方法を見てきました。
重要な点としては、

- 一定距離ずつレイを進める
- オブジェクト内にレイが入ったら、距離関数を元に密度を求める
- オブジェクト内点からライト方向にもサンプリングを行う

というあたりでしょうか。

今回は距離関数を利用して「密度」を計算していましたが、Volume Textureを使ったり、あるいは別の関数を使うなどして形状を他の方法で求めるようにすれば、色々な形の雲（や霧状のオブジェクト）なんかを描くのに使えそうですね。
また、各点におけるサンプリングの手法などは色々なことにも流用できそうなので使いみちがありそうです。

ただ、決して軽い処理ではないので、多様は禁物です。