Unity でリアルな雲を表現するためのシェーダを作成する

# はじめに

Unity で OpenVDB ファイルを読み込んで表示するネイティブプラグインを、オープンソースで開発しています。ソースコードを [github](https://github.com/karasusan/OpenVDBForUnity) で公開しています。
![Introduction.gif](https://qiita-image-store.s3.amazonaws.com/0/74764/a0854860-fc0e-c512-755e-9ff46d999d2a.gif)

## 背景

Oats Studios の [『ADAM』](https://youtu.be/tSDsi2ItktY)、 Unity デモチームの制作した [『Book of the Dead』](https://youtu.be/DDsRfbfnC_A) にみられるように、映像制作に Unity を制作する事例が増えています。 このような状況の背景には、ハードウェア性能の向上はもちろんですが、映像制作に使えるツール群が充実してきたことが理由として挙げられます。
それらのツールの一つに [Alembic Importer](https://github.com/unity3d-jp/AlembicForUnity) がありますが、これは [Alembic フォーマット](http://www.alembic.io)を Unity で利用できるようにすることで、既存のUnityの機能では困難だった表現を可能にしています（[参考](https://blogs.unity3d.com/jp/2017/12/04/adam-the-evolution-of-alembic-support-in-unity)）。映像制作パイプラインで広く利用されているフォーマットや手法を Unity に対応させることで、映像表現の幅を広げた事例といえます。


## OpenVDB
<img src="https://qiita-image-store.s3.amazonaws.com/0/74764/c9836c1a-885b-00f6-c719-2d4b61cbc4a8.jpeg" width=50%>


[OpenVDB](http://www.openvdb.org) は、2012年から主に DreamWorks Animation で開発されている、ボリュームデータを扱うライブラリです。Arnold、VRay、Renderman　などの主要なレンダラでサポートされ、ボリュームデータの業界標準フォーマットとして確立しています。
OpenVDBForUnity はこのフォーマットを Unity でインポートすることで、爆発や雲をリアルタイムレンダリングでリアルに表現することを目指しています。

<!-- 
ボリュームデータの DCC ツールとして SideFX 社の Houdini がデファクト
現在もオープンソースとして開発が続けられています。
-->

## 文章の内容

OpenVDBForUnity の開発で、以下のような知見を得ることができました。

- C++11 でのネイティブプラグイン開発
- [TravisCI](https://travis-ci.org)、[AppVeyor](https://www.appveyor.com/) を用いたオープンソースのCI環境構築
- [CMake](https://cmake.org) 、[Conan](https://conan.io) を用いたマルチプラットフォームビルドパイプライン
- ボリュームレンダリングのシェーダ開発

今回は、これらの中でも特にボリュームレンダリングのシェーダ開発について説明します。

# 解説

OpenVDBForUnity において、ボリュームデータがシェーダに渡るまでの流れを単純に表すと、以下のようになります。

1. OpenVDB ライブラリを Unity のネイティブプラグインとしてインポートする
2. プラグインで vdbファイルからボリュームデータを読み込む
3. ボリュームデータを Texture3D として展開して Unity 側のメモリに転送する
4. ボリュームサイズにあわせたキューブを生成してマテリアルを貼り付ける
4. Texture3D をフラグメントシェーダでレイマーチングしてレンダリングする

今回説明するシェーダ開発は、最後の工程のレンダリングの部分になります。

## 開発環境

開発端末には MacBook Pro を利用していますが、ボリュームレンダリングをリアルタイムに行うには GPU パフォーマンスが要求されるため、開発には eGPU を利用しています。 以下に構成を示します。

- MacBook Pro (15-inch Mid 2017)
- Sonnet eGFX Breakaway Box 650W
- AMD Radeon RX VEGA 64 GPU (GV-RXVEGA64GAMING OC-8GD)
- Unity 2018.2

eGPU Box は市場に出回り始めていますが、 Apple では Sonnet 製品を[推奨](https://support.apple.com/ja-jp/HT208544) しているので、無難に電源容量が650Wのものを利用しています。グラフィックボードは Sonnet 社が公開している [互換性リスト](http://siroanko.hatenablog.com/entry/2018/04/17/003000) を確認してから決めています。
[GeekBench4](https://www.geekbench.com) での計測結果では、Discrete GPU(AMD Radeon Pro 560)との差は3.9倍、内蔵GPU(Intel HD Graphics 630)との差は6.3倍となっています。

![Screen Shot 2018-08-19 at 10.11.03 AM.png](https://qiita-image-store.s3.amazonaws.com/0/74764/eb67c003-3a0f-ced5-1aa4-a659aa681c66.png)

## サンプルデータ

OpenVDBForUnity の開発には VDBファイル（拡張子は **vdb**）が必要になるため、以下のサイトからダウンロードして利用しています。

- Walt Disney Animation Studios (https://www.disneyanimation.com/technology/datasets)
- OpenVDB official site (http://www.openvdb.org/download)

## シェーダファイルの構成

[VolumeStandard.shader](https://github.com/karasusan/OpenVDBForUnity/blob/master/OpenVDBForUnity/Assets/OpenVDB/Shaders/VolumeStandard.shader) で各種シェーダパラメータを定義しています。また、いくつかのプリプロセッサを用意して、表示不具合やパフォーマンスの確認に利用しています。実際の処理は [VolumeStandard.cginc](https://github.com/karasusan/OpenVDBForUnity/blob/master/OpenVDBForUnity/Assets/OpenVDB/Shaders/VolumeStandard.shader) に実装しています。

## シェーダの処理の流れ

フラグメントシェーダの処理は以下のような順序になります。処理順序に沿って解説していきます。

1. ワールド座標系からオブジェクト座標系への変換
2. カメラ位置の内外判定
3. レイの開始位置とサンプリング距離のアラインメント
4. レイの経路に他のオブジェクトがあった場合にレイの終端位置を計算
5. レイマーチング処理
6. レイが何にもぶつからずに通過したときはマスク
7. レイマーチングの結果を格納

ここでは便宜的に、レイマーチングより前の処理を前処理、後の処理を後処理として説明します。

## 前処理

前処理では、レイマーチングに必要なレイの開始位置と終了位置を計算し、１ステップ毎のベクトルを求めます。

### オブジェクト座標系への変換

レイマーチングにはスクリーン座標系とオブジェクト座標系での実装がありますが（[参考記事](http://tips.hecomi.com/entry/2016/09/26/014539)）、このシェーダではオブジェクト座標系で実装しています。頂点シェーダでワールド座標系の位置を計算し、ピクセルシェーダで Localize 関数に変換することでオブジェクト座標系の位置を取得しています。

```csharp:Utils.cginc
inline float3 Localize(float3 p)
{
    return mul(unity_WorldToObject, float4(p, 1)).xyz;
}
```

ボリュームのメッシュ形状は、シェーダ内の計算が簡単という理由でキューブを利用します。オブジェクト座標系を採用することで、Transform の変更を自由に行うことができます。
![TransformVolume.gif](https://qiita-image-store.s3.amazonaws.com/0/74764/2c226034-4e6a-ecf6-f484-7aa19c990587.gif)

### カメラ位置の内外判定

レイの開始位置は、通常はカメラの視線ベクトルとメッシュの面の交点になります。しかし、カメラがメッシュ内部に入ったときは、カメラの near を初期位置に設定する必要があります。カメラ周りの操作は、 hecomi さんの [Camera.cginc](https://github.com/hecomi/uRaymarching/blob/master/Assets/uRaymarching/Shaders/Include/Camera.cginc) のコードをそのまま使っています。

```csharp:VolumeStandard.cginc
    Ray ray;
    ray.origin = Localize(i.world);

    float3 cameraDir = GetCameraDirection(i.pos);
    float3 cameraPos = GetCameraPosition();
    float3 nearCameraPos = cameraPos + (GetCameraNearClip() + 0.01) * cameraDir;
    nearCameraPos = Localize(nearCameraPos);

    if(IsInnerCube(nearCameraPos))
    {
        ray.origin = nearCameraPos;
    }
```

レイの開始位置を調整することで、カメラがボリュームのキューブ内部に入っても正しくレンダリングされるようになります。
![CameraInsideCube.gif](https://qiita-image-store.s3.amazonaws.com/0/74764/50d20490-5151-dc09-e9e4-04fa1ae25b83.gif)

### サンプル位置のアラインメント

`ray.origin` を、そのままレイトレーシングの初期位置として利用してもレンダリングできますが、モアレのような表示が目立ちます（アーティファクトと呼んでいます）。この原因はレイトレーシングのサンプリング位置が近傍とズレていることで生じています。カメラ方向の垂直平面上をサンプリングすることで解決します。

ワールド座標系におけるレイの方向（`cameraDir`）とカメラ方向（`cameraForward`） の内積をとることで余弦（cos） をとります。

```math
\vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}|\cos\theta
```

$\vec{a}$ と $\vec{b}$　がともに単位ベクトルのとき

```math
\vec{a} \cdot \vec{b} = \cos\theta
```

カメラ方向のサンプリングの距離（`_StepDistance`）を基準にすると、 レイ方向のサンプリング距離（`stepDist`）は三角関数で求められます。

```math
\frac{|\vec{a}|}{|\vec{b}|} = \cos\theta
\\
|\vec{b}| = \frac{|\vec{a}|}{\cos\theta}

```

```csharp:VolumeStandard.cginc
    float3 cameraForward = GetCameraForward();
    float stepDistRatio = 1.0 / dot(cameraDir, cameraForward);
    float stepDist = _StepDistance * stepDistRatio;
```


`fmod` 関数でカメラ位置から初期位置までの距離（`cameraDist`） を `stepDist` で割った余りを求めます。

```csharp:VolumeStandard.cginc
    float cameraDist = length(i.world - cameraPos);
    float startOffset = stepDist - fmod(cameraDist, stepDist);

    // sampling parameter (start, end, stepcount)
    float3 start = ray.origin + mul((float3x3) unity_WorldToObject, cameraDir * startOffset);
```

### 他オブジェクトとの交差判定

レイマーチングでは、他のオブジェクトがボリュームに含まれているときの交差判定を記述する必要があります。交差判定のあり/なしを`ENABLE_TRACE_DISTANCE_LIMITED` プリプロセッサの有無で比較することができます。
![TraceDistanceLimited.gif](https://qiita-image-store.s3.amazonaws.com/0/74764/a75d8dca-2419-588c-148f-79354b932153.gif)
他オブジェクトとの交差判定を行うために深度テクスチャ（`_CameraDepthTexture`）を使用します。

```csharp:VolumeStandard.cginc
    float2 uv = i.screenPos.xy / i.screenPos.w;
    float sceneDepth = LinearEyeDepth(SAMPLE_DEPTH_TEXTURE(_CameraDepthTexture, uv));

    float tfar2 = length(ray.origin - Localize(sceneDepth * cameraDir + cameraPos));
    end = ray.origin + ray.dir * min(tfar, tfar2);
```

ここで `i.screenPos` は、頂点シェーダ側での `ComputeScreenPos` の計算結果です。シーンの深度からレイの到達するまでの長さを計算し、`min(tfar, tfar2)` の比較をすることで、距離が近い方をレイの終了位置までの距離として利用します。

### ステップ単位ベクトルの計算

レイの開始位置と終了位置と１ステップ毎の距離が分かれば、１ステップ毎のベクトル(`ds`) とステップ数（`stepCount`）が求められます。

```csharp:VolumeStandard.cginc
    float dist = length(end - start);
    half stepCount = dist / stepDist;
    float3 ds = ray.dir * stepDist;
```

## レイマーチング

これまで説明してきた処理は、主にレイの開始位置と終了位置を求めるための計算で、実際のシェーディングに関する処理はこのレイマーチングの処理で行われます。
レイマーチングは、レイの開始位置から終了位置まで一定の間隔を刻みながら、ボクセルの値をサンプリングします。ここのロジックは繰り返しの処理になるので、処理が複雑になるとパフォーマンスに大きく影響します。

### サンプリング

現在位置(`p`)からボクセルをサンプリングします。サンプル結果が `0.01` 未満の場合は残りの処理をスキップします。`p += ds` でレイのステップを進めます。ステップ数をオーバーするか、あるいは透過率（`transmittance`）が `0.01` を越えた場合にループを抜けます。

```csharp:VolumeStandard.cginc
    for (int iter = 0; iter < ITERATIONS; iter++)
    {
        float3 uv = GetUV(p);
        float cursample = SampleVolume(uv);

        if(cursample > 0.01)
        {
            ...
        }
        p += ds;

        if(iter >= stepCount)
        {
            break;
        }

        if (transmittance < 0.01)
        {
            break;
        }
    }
```

サンプル結果が 0.01 以上だった場合、以下の処理を行います。

- ライト方向へのサンプリング
- 環境光のシミュレートのサンプリング

### ライト方向へのサンプリング

ディレクショナルライト（１灯のみ）の影響を受けるように実装します。ライト方向にサンプリングすることで、そのボクセルにおける影の距離を計算します。透過率のしきい値を越えるか、ボックスの領域を越えたときに計算を早期終了します。

```csharp:VolumeStandard.cginc

            float3 lpos = p;
            float shadowdist = 0;                        

            for (int s = 0; s < _ShadowSteps; s++)
            {
                lpos += lightVec;
                float3 luv = GetUV(lpos);
                float lsample = SampleVolume(saturate(luv));

                shadowdist += lsample;

                float3 shadowboxtest = floor( 0.5 + ( abs( 0.5 - luv ) ) );
                float exitshadowbox = shadowboxtest .x + shadowboxtest .y + shadowboxtest .z;

                // check to exit shadow box
                if(shadowdist > shadowthreshold || exitshadowbox >= 1)
                {
                    break;
                }
            }
```

`ENABLE_DIRECTIONAL_LIGHT` のオン/オフ でライトの影響あり/なしを比較できます。
![Screen Shot 2018-08-19 at 7.44.12 PM.png](https://qiita-image-store.s3.amazonaws.com/0/74764/b5ba9838-3629-06c1-642d-187da5129b70.png)

### 不透明度の計算

光学における[透過率](https://en.wikipedia.org/wiki/Transmittance)の公式、ランベルト・ベールの法則を用いて透過率を計算しています。媒質の吸収係数（ここでは密度）と経路長から透過率が求まります。`shadowterm` は光源方向に対する透過率で、透過率が高いほど明るく（光が届く）ようになります。

```csharp:VolumeStandard.cginc
            curdensity = saturate(cursample * _Intensity);
            float3 shadowterm = exp(-shadowdist * shadowDensity);
            float3 absorbedlight = shadowterm * curdensity;
            lightenergy += absorbedlight * transmittance;
            transmittance *= 1-curdensity;
```

### 環境光

環境光をシミュレートするために、オフセットサンプルを追加することで実現します。

```csharp:VolumeStandard.cginc
            shadowdist = 0;

            float3 luv = uv + float3(0,0,0.05);
            shadowdist += SampleVolume(saturate(luv));
            luv = uv + float3(0,0,0.1);
            shadowdist += SampleVolume(saturate(luv));
            luv = uv + float3(0,0,0.2);
            shadowdist += SampleVolume(saturate(luv));
            lightenergy += exp(-shadowdist * _AmbientDensity) * curdensity * _AmbientColor * transmittance;
```

`ENABLE_AMBIENT_LIGHT` プリプロセッサマクロの有無で、環境光の影響の比較ができます。
![Screen Shot 2018-08-19 at 11.58.24 PM.png](https://qiita-image-store.s3.amazonaws.com/0/74764/c06c9f4b-0375-de94-30ba-3a9ca7370195.png)

## 後処理

レイマーチングで得られた計算結果を、フラグメントシェーダの出力として整形します。

### 不要なピクセルをマスク

一度もサンプリング処理を行わなかったピクセルは、不要な描画を避けるためにマスクします。`clip(-1)` を呼び出せばマスクされます。

```csharp:VolumeStandard.cginc
    if(depthtest)
    {
        clip(-1);
    }
```

### レイマーチングの結果を格納

フラグメントシェーダの戻り値の構造体は以下のように定義しています。

```csharp:VolumeStandard.cginc
struct fragOutput
{
    fixed4 color : SV_Target0;
    float depth : SV_Depth;
};
```

深度の計算については ComputeDepth 関数を [Util.cginc](https://github.com/karasusan/OpenVDBForUnity/blob/master/OpenVDBForUnity/Assets/OpenVDB/Shaders/Utils.cginc) に実装しています（[参考](http://i-saint.hatenablog.com/entry/2015/03/18/194818)）。

```csharp:VolumeStandard.cginc
    fragOutput o;
    o.color = float4(lightenergy, 1-transmittance);
    o.depth = ComputeDepth(UnityObjectToClipPos(float4(depth, 1.0)));
    return o;
```

# おわりに

ボリュームレンダリングのシェーダについて大まかに解説しました。今後対応を予定している項目を列挙しています。また、ボリュームレンダリングについて、より詳しい内容を学習したいという方のために、参考資料へのリンクを列挙しています。

## 今後の対応

- density 以外のパラメータへの対応 （velocity, heat, temperature）
- ShaderGraph への対応
- Timeline によるアニメーション制御
- Houdini を利用した作例の紹介

## 参考資料

以下にシェーダ開発で参考にした資料を列挙します。特にボリュームレンダリングについては、多くのアイデアを [Shader Bits](https://shaderbits.com/blog)の[レイマーチングの記事](http://shaderbits.com/blog/creating-volumetric-ray-marcher) から得ています。また、Unity でのレイマーチングの実装については、凹みさんの[ブログ](http://tips.hecomi.com)が日本語の記事でよくまとめられています。

- [Creating a Volumetric Ray Marcher - Shader Bits] (https://shaderbits.com/blog/creating-volumetric-ray-marcher)
- [Unity でボリュームレンダリングをしてみる - 凹みtips] (http://tips.hecomi.com/entry/2018/01/08/211937)
- [object space raymarching - primitive: blog](http://i-saint.hatenablog.com/entry/2015/08/24/225254)
- [Raymarching Distance Fields: Concepts and Implementation in Unity - adrian's soapbox](http://flafla2.github.io/2016/10/01/raymarching.html)