GLSLでボリュームレイマーチング

GLSLでボリュームレイマーチングを実装してみました。理解するまでだいぶ悩んだので備忘録として記事にしておきます。

今回実装したサンプルは以下から見ることができます。
https://neort.io/art/bj8438k3p9f9psc9p4d0
![result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/0ab2e868-fa99-0d13-7198-8790e001f3ca.png)

以下のマクロで密度関数のベースとなる形状を変更することができます。

```glsl
// 0: box, 1: sphere: 2: torus
#define DENSITY_BASIC_SHAPE 1
```


## アルゴリズム概説


今回実装したボリュームレイマーチングの実装は以下の記事を主に参考にしています。
[Creating a Volumetric Ray Marcher](https://shaderbits.com/blog/creating-volumetric-ray-marcher/)

アルゴリズムを疑似コードっぽく書くと以下のようになります。

```
レイの開始位置を求める
色 = vec3(0.0)
透明度 = 1.0
for (レイが終了位置に達していない)　{
  レイの現在位置の密度をサンプルする
  レイの現在位置に届く光をサンプルする
  色 += 光 * 密度 * 透明度
  透明度 *= 1.0 - 密度
  if (透明度がほぼ0) {
    break;
  }
  レイを一定距離だけ進める
}
不透明度 = 1.0 - 透明度
```

ボリュームレイマーチングでは、カメラから出たレイを一定距離ずつ進めていきます。レイを進めた各地点では、密度(density)とそこに届く光をサンプルして色に加算しつつ、透明度を減らしていきます。透明度がほぼ0になるとその先にあるボリュームは見えないのでレイマーチングを終了します。レイが終了位置に達した場合にもレイマーチングを終了します。

`色 += 光 * 密度 * 透明度`は一見複雑ですが、各要素ごとに見ていくとわかりやすいです。まず、その地点に届く光が大きいほど色に影響します。もし光が届かない(光=0)ならば、その地点は色に影響しません。また、密度が大きいほど多くの光が拡散されるため、色に影響します。もし密度が0ならばその地点で光は拡散されません。最後に透明度が大きいほど、その地点がカメラからよく見えるので影響が大きくなります。透明度が0ならばその地点はカメラから見えないので影響しません。

## バウンディングボックス

ボリュームレイマーチングでは、ボリュームのある箇所のみをサンプリングできると効率がよくなります。今回はボリュームのある箇所を軸平行境界ボックス(Axis-Aligned Bounding Box; AABB)で限定しています。AABBは面がXYZ軸にそれぞれ平行な直方体です。AABBとレイの交点からレイマーチングを開始して、レイがAABBを抜けたらレイマーチングを終了するようにしています。



```glsl
#define RAYMARCH_ITERATIONS 48
#define BOUNDING_BOX_SIZE 6.0
#define BOUNDING_BOX_OFFSET vec3(0.0)

vec4 raymarch(vec3 rayOrigin, vec3 rayDirection, float tmin, float tmax) {
    // 今回のAABBは立方体なのでレイマーチングの最大距離は対角線の長さになる
    // 最大距離をレイマーチング数で割ることで、1ステップあたりのレイを進める距離を求める
    float raymarchSize = sqrt(3.0) * BOUNDING_BOX_SIZE / float(RAYMARCH_ITERATIONS);
    // レイがAABBを脱出するまでのイテレーション数を求める
    int maxRaymarchIteration = int((tmax - tmin) / raymarchSize) + 1;
    // レイの開始点はAABBとの交点のうち近いほうになる
    vec3 position = rayOrigin + tmin * rayDirection - BOUNDING_BOX_OFFSET;
    ...
    for (int ri = 0; ri < RAYMARCH_ITERATIONS; ri++) {
        if (ri >= maxRaymarchIteration) { // レイがAABBを脱出したら途中終了する
            break;
        }
        ...
    }
    ...
}

// AABBにヒットしたときtrueを返し、ヒットしないときfalseを返す
// ヒットしたときレイマーチの開始点がro + tmin * rd、終了点がro + tmax * rdになる
bool aabb(vec3 ro, vec3 rd, vec3 corner0, vec3 corner1, inout float tmin, inout float tmax) {
    for (int i = 0; i < 3; i++) {
        float t0 = (corner0[i] - ro[i]) / rd[i];
        float t1 = (corner1[i] - ro[i]) / rd[i];
        tmin = max(tmin, min(t0, t1));
        tmax = min(tmax, max(t0, t1));
        if (tmax <= tmin) {
            return false;
        }
    }
    return true;
}

void main(void) {
    ...
    // AABBの2つの角を定義する
    vec3 corner0 = vec3(-0.5 * BOUNDING_BOX_SIZE + BOUNDING_BOX_OFFSET);
    vec3 corner1 = vec3(0.5 * BOUNDING_BOX_SIZE + BOUNDING_BOX_OFFSET);
    float tmin = 0.0;
    float tmax = 1e6;
    if (aabb(ro, rd, corner0, corner1, tmin, tmax)) { // AABBにヒットしたらレイマーチングを行う
        vec4 res = raymarch(ro, rd, tmin, tmax);
        ...
    }
    ...
}

```

雲海のように画面全体に雲が存在する表現ではカメラ位置からサンプリングを開始すればいいので、バウンディングボックスは不要です。

## 不透明度

まずは、ライティングを考慮せずに不透明度のみを計算したいと思います。

```glsl
vec4 raymarch(vec3 rayOrigin, vec3 rayDirection, float tmin, float tmax) {
    float raymarchSize = sqrt(2.0) * BOUNDING_BOX_SIZE / float(RAYMARCH_ITERATIONS);
    float densityScale = DENSITY_INTENSITY * raymarchSize;
    vec3 rayStep = rayDirection * raymarchSize;
    vec3 color = vec3(0.0); // ここでは使用しない
    float transmittance = 1.0; // 透明度
    for (int ri = 0; ri < RAYMARCH_ITERATIONS; ri++) {
        ...
        float density = sampleDensity(position);
        if (density > 0.001) { // 密度が一定以下の場合、効率化するために計算を省略する
            density = saturate(density * densityScale);
            transmittance *= 1.0 - density;
        }
        if (transmittance < 0.001) { // 透明度が0に近づいたら終了する
            break;
        }
        position += rayStep; // レイを進める
    }
    return vec4(color, 1.0 - transmittance);
}


void main(void) {
    ...
    vec3 color = background(rd);
    ...
    if (aabb(ro, rd, corner0, corner1, tmin, tmax)) {
        vec4 res = raymarch(ro, rd, tmin, tmax);
        // 透明度が残っている場合は背景が透けるようにする
        color = res.xyz + (1.0 - res.w) * color;
    }
    ...
}
```

サンプルでは以下のようにマクロをコメントアウトすると、不透明度のみを考慮してレンダリングされます。

```glsl
// #define WITH_DIRECTIONAL_LIGHT
// #define WITH_AMBIENT_LIGTHT
```

![opacity.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/541210c3-38cc-00b6-5c0d-6d1eff355125.png)


## 平行光源によるライティング

次に平行光源によるライティングを考慮したいと思います。各サンプル点でさらに平行光源方向にレイマーチングを行うことで光方向のボリュームの厚さを調べ、それにより光がどれぐらい吸収されるかを求めます。吸収されなかった光の拡散を色に追加します。

光の吸収は以下の式であらわされる[ランベルト・ベールの法則](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3%E3%83%99%E3%83%AB%E3%83%88%E3%83%BB%E3%83%99%E3%83%BC%E3%83%AB%E3%81%AE%E6%B3%95%E5%89%87)をもとに行います。$A$が吸収度、$t$が光の移動距離、$d$が密度を表しています。

```math
A = e^{-td}
```

レイマーチングでの1ステップあたりのレイの移動距離を$t'$とすると以下のように近似できます。

```math
A = e^{-td} \approx e^{-\sum_{i=0}^{N} t_{i}d_{i}} = e^{-t'\sum_{i=0}^{N} d_{i}}
```

```glsl
#define SHADOW_ITERATIONS 8
#define DENSITY_INTENSITY 2.0

const vec3 DENSITY_COLOR = vec3(0.8, 0.9, 1.0);
const vec3 ABSORPTION_INTENSITY = vec3(0.5, 0.8, 0.7) * 0.5;
const vec3 DIRECTIONAL_LIGHT_DIR = normalize(vec3(0.5, 1.0, 0.5));
const vec3 DIRECTIONAL_LIGHT_COLOR = vec3(1.0, 1.0, 0.8) * 1.0;

vec4 raymarch(vec3 rayOrigin, vec3 rayDirection, float tmin, float tmax) {
    float raymarchSize = sqrt(3.0) * BOUNDING_BOX_SIZE / float(RAYMARCH_ITERATIONS);
    float shadowSize = SHADOW_LENGTH / float(SHADOW_ITERATIONS);
    float densityScale = DENSITY_INTENSITY * raymarchSize;
    vec3 shadowScale = ABSORPTION_INTENSITY * shadowSize;
    int maxRaymarchIteration = int((tmax - tmin) / raymarchSize) + 1;
    vec3 position = rayOrigin + tmin * rayDirection - BOUNDING_BOX_OFFSET;
    vec3 rayStep = rayDirection * raymarchSize;
    vec3 shadowStep = DIRECTIONAL_LIGHT_DIR * shadowSize;
    vec3 color = vec3(0.0);
    float transmittance = 1.0;
    for (int ri = 0; ri < RAYMARCH_ITERATIONS; ri++) {
        if (ri >= maxRaymarchIteration) {
            break;
        }
        float density = sampleDensity(position);
        if (density > 0.001) {
            density = saturate(density * densityScale);

            vec3 shadowPosition = position;
            float shadowDensity = 0.0;
            // 平行光源方向にレイマーチングをおこない、吸収による光の減衰を調べる
            for (int si = 0; si < SHADOW_ITERATIONS; si++) {
                shadowPosition += shadowStep;
                shadowDensity +=  sampleDensity(shadowPosition);
            }
            vec3 attenuation = exp(-shadowDensity * shadowScale);
            vec3 attenuatedLight = DIRECTIONAL_LIGHT_COLOR * attenuation; // 吸収により減衰した光
            color += DENSITY_COLOR * attenuatedLight * transmittance * density; // 吸収されなかった光は拡散する
            ...
            transmittance *= 1.0 - density;
        }
        if (transmittance < 0.001) 
           break;
        }
        position += rayStep;
    }
    return vec4(color, 1.0 - transmittance);
}
```

マクロを以下のようにすると平行光源のみを考慮したレンダリングを行います。

```
#define WITH_DIRECTIONAL_LIGHT
// #define WITH_AMBIENT_LIGTHT
```

![directional.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/a1372a4a-60b5-969f-47ce-e32e7b3b2f01.png)


## 環境光によるライティング

最後に環境光を考慮したいと思います。適当な方向にレイを投げて、それを環境光として加算します。環境光は平行光源のシャドウ部分に陰影を出すために使うので、環境光とは反対方向にレイを投げるといいと思います。

```glsl
const vec3 AMBIENT_LIGHT_DIR = normalize(vec3(0.0, -1.0, 0.0));
const vec3 AMBIENT_LIGHT_COLOR = vec3(0.5, 0.7, 1.0) * 0.2;

vec4 raymarch(vec3 rayOrigin, vec3 rayDirection, float tmin, float tmax) {
    ...
    float transmittance = 1.0;
    for (int ri = 0; ri < RAYMARCH_ITERATIONS; ri++) {
        if (ri >= maxRaymarchIteration) {
            break;
        }
        float density = sampleDensity(position);
        if (density > 0.001) {
            density = saturate(density * densityScale);
            ...
            { // 環境光の計算（計算自体は平行光源とほぼ同じ）
                float shadowDensity = 0.0;
                vec3 shadowPosition = position + AMBIENT_LIGHT_DIR * 0.05;
                shadowDensity += sampleDensity(shadowPosition) * 0.05;
                shadowPosition = position + AMBIENT_LIGHT_DIR * 0.1;
                shadowDensity += sampleDensity(shadowPosition) * 0.05;
                shadowPosition = position + AMBIENT_LIGHT_DIR * 0.2;
                shadowDensity += sampleDensity(shadowPosition) * 0.1;
                float attenuation = exp(-shadowDensity * AMBIENT_INTENSITY); // attenuated by absorption
                vec3 attenuatedLight = AMBIENT_LIGHT_COLOR * attenuation;
                color += DENSITY_COLOR * attenuatedLight * transmittance * density; // out-scattering
            }
            transmittance *= 1.0 - density;
        }
        ...
    }
    return vec4(color, 1.0 - transmittance);
}
```

環境光のみでレンダリングする場合は、サンプルのマクロを以下のようにしてください。

```glsl
// #define WITH_DIRECTIONAL_LIGHT
#define WITH_AMBIENT_LIGTHT
```

![ambient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/21b73c8f-1551-6dcd-9fb8-59d5b4c52aa7.png)

## 結果

平行光源と環境を両方組み合わせると以下のようになります(最初の画像と同じです)。
この画像だとわかりづらいですが、平行光源の影になっている部分を見ると環境光の有無の影響がわかると思います。
![result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/0ab2e868-fa99-0d13-7198-8790e001f3ca.png)


以上、GLSLでボリュームレイマーチングを行う方法について解説しました。私自身の理解が足りていない部分もあるので、間違っている箇所があれば指摘していただけると幸いです。

## 参考

* [Creating a Volumetric Ray Marcher](https://shaderbits.com/blog/creating-volumetric-ray-marcher/)
* [Unity でリアルな雲を表現するためのシェーダを作成する \- Qiita](https://qiita.com/karasusan/items/ce8891f5c4b63d7d06ba)
* [Cloudy Shapes](https://www.shadertoy.com/view/WdXGRj)