2019年12月24日火曜日

clspvを用いてopenclカーネルをvulkanで実行する

Clspvを用いてOpenCLカーネルをVulkanで実行する

こんにちは。このページは レイトレ Advent Calendar 2019 の19日目の記事です。この記事では以下の事について紹介します、

  • Vulkan SDK と Clspv のインストール
  • Clspv で OpenCL カーネルをビルド
  • ビルドした カーネルを Vulkan で実行する際の注意事項

この記事では Vulkan 仕様や API の使い方については説明しません。Vulkan Compute Shader の使い方は Clspv を使った場合でも特別な事は無いため、他の資料などを参考にしてもらえればと思います。

Clspv を使った Vulkan のサンプルコードは github にアップロードしていますので参考にしてください。

Vulkan SDK のインストール

まずは、下記のページから最新の Vulkan SDK をダウンロードします。

SDKのインストール方法は Windows、Linux、macOS で異なります。

Windows

ダウンロードしたインストーラーを起動しそのままインストールします。インストールが完了したら VULKAN_SDK という環境変数が作られ、Vulkan SDK へのパスが設定されます。

Vulkan environment variable

Linux

ダウンロードした tar.gz を展開します。展開したディレクトリ内に setup-env.sh というファイルがあり、これをシェルに反映させると、SDK内のツールやライブラリへのパスを自動的に設定してくれます。また、 VULKAN_SDK という環境変数がつくられて、Vulkan SDKへのパスが設定されます。

% source  ${path-to-vulkan-sdk}/setup-env.sh
% echo  ${VULKAN_SDK}
${path-to-vulkan-sdk}/x86_64

macOS

ダウンロードした tar.gz を展開します。Linux版と違い展開したディレクトリ内に setup-env.sh が無いため以下のようなスクリプトを用意します。

export VULKAN_SDK="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &&  pwd )"/macOS
export PATH="${VULKAN_SDK}/bin:${PATH}"
export DYLD_LIBRARY_PATH="${VULKAN_SDK}/lib:${DYLD_LIBRARY_PATH}"
export VK_LAYER_PATH="${VULKAN_SDK}/etc/vulkan/explicit_layer.d"
export VK_ICD_FILENAMES="${VULKAN_SDK}/etc/vulkan/icd.d/MoltenVK_icd.json"

このスクリプトを展開したディレクトリ内に置きシェルに反映させると、SDK内のツールやライブラリへのパスを自動的に設定してくれます。

Vulkan SDK の動作確認

Vulkan SDK をインストールしたのできちんと動作するか確認します。動作確認には vkvia を用います。vkvia は簡単なVulkanアプリケーションをいくつか動かして、その時のVulkanの状態をキャプチャーします。

vkvia1

vkvia${VULKAN_SDK}/bin/ にあります。実行すると実行ディレクトリに vkvia.html が作成され、ブラウザで開くとキャプチャーしたVulkanの状態を見ることができます。例えば、 Physical Devices を開くとVulkanで使用できるハードウェアデバイスの情報を見ることができます。

vkvia2

Clspv のインストール

Clspvにはリリースバージョンの実行ファイルやソースコードはまだありません。github の google/clspv から master ブランチを pull してビルドします。ビルドするのに必要なツールはREADMEに書いてあるように、

  • CMake
  • Python
  • C++ compiler
  • git

になります。まずはClspvが依存する外部プロジェクトのソースコードをダウンロードします。

% python utils/fetch_sources.py

ダウンロードした後は CMake を実行して Clspv のビルドの準備をします。Windows の場合は (Visual Studio 2019 を使う場合)、

% cmake -G"Visual Studio 16 2019" -Ax64 -Thost=x64 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${path-to-dir} ..

${path-to-dir} にはインストール先のパスが入ります。
Linux や macOS の場合は、

% cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${path-to-install} ..

です。 CMake でジェネレーターファイルを作成したら Clspv をビルドします。

% cmake --build . --config Release --target install --parallel 8

ビルドが完了すると ${path-to-dir}clspv という実行ファイルがインストールされます。

Clspv で OpenCL コードをビルド

Clspv は OpenCL C 1.2 のコードをビルドすることができますが、Vulkan が扱うことができない機能については Clspv 側でビルドができないように制限がかけられています。制限内容の詳細は OpenCL C Restrictions を読んで下さい。いくつかピックアップすると、

  • 8や16長のベクトルは使用できない ( int8 など)
  • double は使用できない
  • cbrtpownrintsincos など使用できない Math 関数がある
  • printf は使用できない

などがあります。ただ、これらの制限は将来 Vulkan の仕様のアップデートに合わせて緩和されるかもしれません。

では、実際に Clspv で以下の OpenCL C コードをビルドしてみます。

// vulkan_clspv_test0.cl
__kernel void test(__global const float4* inputs, __global float4* outputs)
{
  const uint index = get_global_id(0);
  outputs[index] = sqrt(inputs[index]);
}

コマンドは以下のようになります。

% clspv -o=vulkan_clspv_test0.spv vulkan_clspv_test0.cl

-o オプションで出力する SPIR-V ファイルのファイル名を指定できます。その他にも最適化オプション -O=3 などがあり、オプション一覧は clspv --help で見ることができます。

OpenCL C++ のビルド

Clspv は実験的ですが OpenCL C++ をサポートします。現状 OpenCL C++ の仕様を完全にサポートしているわけではく、いくつかの違いはあります。例えば、

  • OpenCL C++ Standard Library は使用できない
    • Mathライブラリなどは OpenCL C の組み込み関数スタイルの物を使用する
  • C++14 ではなく C++17 ベース
    • constexpr iffold expression など C++17 の機能も使うことができる

等の違いがあります。Clspv で OpenCL C++ のコードをビルドする場合は --c++ --inline-entry-points オプションを付けてビルドします。例として、以下の OpenCL C++ コードをビルドしてみます。

// vulkan_clspv_test1.cl
namespace clspvtest {

// Forward declaration
float sum(const float x);
float sum(const float4 x);
template <typename... Types> float sum(const Types... args);
class Matrix4x4;
Matrix4x4 operator+(const Matrix4x4& lhs, const Matrix4x4& rhs);

float sum(const float x)
{
    return x;
}

float sum(const float4 x)
{
    const float s = x.x + x.y + x.z + x.w;
    return s;
}

template <typename... Types>
float sum(const Types... args)
{
    // Clspv can build C++17 features
    if constexpr (0 < sizeof...(Types)) {
        const auto s = (args + ...);
        return sum(s);
    }
    else {
        return 0.0f;
    }
}

class Matrix4x4
{
  public:
    Matrix4x4() {}
  
    Matrix4x4(const float4 r1, const float4 r2, const float4 r3, const float4 r4) :
        r1_{r1}, r2_{r2}, r3_{r3}, r4_{r4}
    {
    }

    //! Return the sum of the elements
    float sum() const
    {
        const float s = clspvtest::sum(r1_, r2_, r3_, r4_);
        return s;
    }

    float4 r1_ = float4{0.0f, 0.0f, 0.0f, 0.0f},
           r2_ = float4{0.0f, 0.0f, 0.0f, 0.0f},
           r3_ = float4{0.0f, 0.0f, 0.0f, 0.0f},
           r4_ = float4{0.0f, 0.0f, 0.0f, 0.0f};
};

//! Add corresponding elements of the given two matrices
Matrix4x4 operator+(const Matrix4x4& lhs, const Matrix4x4& rhs)
{
    const Matrix4x4 result{lhs.r1_ + rhs.r1_, lhs.r2_ + rhs.r2_, lhs.r3_ + rhs.r3_, lhs.r4_ + rhs.r4_};
    return result;
}  

} // namespace clspvtest 

__kernel void testSummation(__global clspvtest::Matrix4x4* value, __global float* outputs)
{
    const uint index = get_global_id(0);
    if (0 < index)
        return;

    // Check clspvtest::Matrix4x4
    static_assert(sizeof(clspvtest::Matrix4x4) == 64, "The size of clspvtest::Matrix4x4 isn't 64 bytes.");
    static_assert(alignof(clspvtest::Matrix4x4) == 16, "The alignment of clspvtest::Matrix4x4 isn't 16 bytes.");

    // Constant value
    __constant constexpr float k = 4.0f; // constant value must be in the outermost scope of a kernel
    static_assert(k == 4.0f);  

    outputs[0] = clspvtest::sum();
    outputs[1] = clspvtest::sum(sqrt(k)); // Use OpenCL C style built-in function
    outputs[2] = clspvtest::sum(1.0f, 2.0f, 3.0f);

    // Matrix4x4 in registers
    {
        __private const clspvtest::Matrix4x4 m;
        outputs[3] = m.sum();
    }
    {
        __private const clspvtest::Matrix4x4 m{float4{k, k, k, k},
                                               float4{k, k, k, k},
                                               float4{k, k, k, k},
                                               float4{k, k, k, k}};
        outputs[4] = m.sum();
    }

    // Matrix4x4 in global memory
    {
        __private const clspvtest::Matrix4x4 m1{float4{k, k, k, k},
                                                float4{k, k, k, k},
                                                float4{k, k, k, k},
                                                float4{k, k, k, k}};
        __private const clspvtest::Matrix4x4 m2{float4{0.0f, 1.0f, 2.0f, 3.0f},
                                                float4{4.0f, 5.0f, 6.0f, 7.0f},
                                                float4{8.0f, 9.0f, 10.0f, 11.0f},
                                                float4{12.0f, 13.0f, 14.0f, 15.0f}};
        value[0] = m1 + m2; // Initialize global Matrix4x4
        outputs[5] = value->sum();
    }

    // Matrix4x4 in local memory
    __local clspvtest::Matrix4x4 storage[1]; // local variable must be in the outermost scope of a kernel
    {
        __private const clspvtest::Matrix4x4 m = value[0];
        storage[0] = m;
        outputs[6] = storage->sum();
    }
}

以下のコマンドでビルドします。

% clspv --c++ --inline-entry-points -o=vulkan_clspv_test1.spv vulkan_clspv_test1.cl

ビルドが成功して SPIR-V ファイルが生成されます。実行結果も以下のように、

output[0] = 0
output[1] = 2
output[2] = 6
output[3] = 0
output[4] = 64
output[5] = 184
output[6] = 184

正しく実行されています。( github のサンプルコードをビルドすると実行できます。)

Vulkan で OpenCL カーネルを実行する

Clspv で生成した SPIR-V コードは Vulkan Compute Shader で実行できます。SPIR-V コードの実行方法は通常の方法と変わらないためここでは省略します。Clspv の SPIR-V を実行する上で気をつける点がいくつかあり、詳細については OpenCL C 1.2 Language on Vulkan を参照して下さい。ここでは大事なものをいくつかピックアップして説明します。

DescriptorSet

Clspv でビルドした OpenCL カーネルの引数は、 Vulkan では DescriptorSet にバインドしたバッファーに対応します。基本的に OpenCL カーネルの引数の順番と DescriptorSet にバインドしたバッファーの順番は対応しています。カーネル引数とバッファーの具体的な対応関係を確認したい場合は、Clspv で Descriptor map を出力できます。Descriptor map は Clspv でビルドする際に --descriptorma=<filename>を付けることで出力できます。例えば、 前節の vulkan_clspv_test1.cl の場合は、

% clspv --c++ --inline-entry-points --descriptormap=vulkan_clspv_test1.csv -o=vulkan_clspv_test1.spv vulkan_clspv_test1.cl

とすることで vulkan_clspv_test1.csv というファイル名のCSV形式の Descriptor map が出力されます。CSVファイルの内容は、

kernel,testSummation,arg,value,argOrdinal,0,descriptorSet,0,binding,0,offset,0,argKind,buffer
kernel,testSummation,arg,outputs,argOrdinal,1,descriptorSet,0,binding,1,offset,0,argKind,buffer

のようになります。各行はそれぞれひとつの引数についての情報が載っています。行の内容は、フィールド名と値が交互に記述されています。行の内容は、

  • kernel: OpenCLカーネル名
  • arg: 引数名
  • argOrdinal: 引数のポジション
  • descriptorSet: DescriptorSet のインデックス
  • binding: 使用するバインドされたバッファーのポジション
  • offset: バッファーにアクセスする際のアドレスのオフセット
  • argKind: バッファの種類

となります。この、binding の部分を見ればカーネルの引数がどのバッファーに対応するかわかります。

Work-Group Size

OpenCL では Work-Group のサイズは clEnqueueNDRangeKernel を呼ぶ時に引数として指定していました。Clspv でビルドしたカーネルの Work-Group サイズを指定する場合は、reqd_work_group_size を使ってカーネルの属性として指定するか、Specialization Constants を設定して指定する方法があります。Vulkan を動かす GPU が固定な場合は reqd_work_group_size を使用してハードコードする方法で良いですが、動かすGPUが複数あり、GPUによって Work-Group サイズを変える場合は Specialization constants で指定する方法が良いです。Specialization Constants は Compute Shader のパイプラインを作成する時に設定します。

std::array<uint, 3> work_group_size{{x, y, z}}; // X、Y、Z次元の Work-Group サイズ
// X、Y、Zの Work-Group サイズを設定するため、
// それぞれの specialization constant の ID やデータサイズを設定する
std::array<vk::SpecializationMapEntry, 3> entries;
for (std::size_t i = 0; i < entries.size(); ++i) {
    entries[i].constantID = static_cast<uint>(i);
    entries[i].offset = static_cast<uint>(i * sizeof(uint));
    entries[i].size = sizeof(uint);
}
const  vk::SpecializationInfo info{3, entries.data(), 3 * sizeof(uint), work_group_size.data()};
// Work-Group サイズの情報を含めて パイプラインを作成する
vk::ShaderModule module = /* Clspv でビルドした SPIR-V を読み込んだモジュール */
const  char* kernel_name = /* OpenCL カーネル名 */
const  vk::PipelineShaderStageCreateInfo shader_stage_create_info{
    vk::PipelineShaderStageCreateFlags{},
    vk::ShaderStageFlagBits::eCompute,
    module,
    kernel_name,
    &info};
vk::PipelineLayout pipeline_layout = /* 作成した PipelineLayout */
const  vk::ComputePipelineCreateInfo create_info{
    vk::PipelineCreateFlags{},
    shader_stage_create_info,
    pipeline_layout};
vk::Device device = /* 作成した Device */
vk::Pipeline pipeline = device.createComputePipelines(vk::PipelineCache{}, create_info);

例えば、Work-Group の合計サイズを 64 で設定する場合、1次元の Work-Group にする場合は、X=64Y=1Z=1 (X * Y * Z = 64) と指定すると、OpenCL 上では、

get_local_size(0); // X = 64
get_local_size(1); // Y = 1
get_local_size(2); // Z = 1

と反映されます。2次元の Work-Group にする場合は、X=8Y=8Z=1 (X * Y * Z = 64) と指定すると、OpenCL 上では、

get_local_size(0); // X = 8
get_local_size(1); // Y = 8
get_local_size(2); // Z = 1

と反映されます。

おわりに

この記事では Vulkan SDK と Clspv の導入を簡単に説明しました。Clspv の登場によって OpenCL コードを Vulkan で動かすことが可能になりました。OpenCL は C言語に慣れている方であれば書きやすい言語だと思いますが、OpenCL ドライバのバージョンが 1.2 で止まっていたりと実行環境周りに不安を覚えている方もいると思います。そのような場合は、 Clspv を使って Vulkan を始めてみるのも面白いと思います。

Clspv はまだ開発中のプロトタイプであるため色々と不具合が見つかるかもしれません。その時は github で issue を作って報告してあげて下さい。

2018年7月8日日曜日

Blenderから自作レンダラー向けに設定をエクスポート

blender_exporting_settings

Blenderから自作レンダラー向けに設定をエクスポート

こんにちは。 このページは レイトレ合宿6 の記事です.

今年もレイトレ合宿の日が近くなってきて どんなシーンをレンダリングしようかと考えている頃ではないでしょうか。 そこで、 1 Blenderで扱うシーンの設定を自作のレンダラー向けにエクスポートできるように してみます。 Blenderからエクスポートすることの利点は、

  • Blenderという強力なGUIツールを用いて効率良くシーンを作成できる
  • Cyclesレンダラーなどの他のエンジンと同じ設定でレンダリングすることで、レンダリング結果を比較できる

などが挙げられます。

デモファイルのダウンロード

設定をエクスポートするに当たって、 Blenderの公式サイトで配布されている ヘリコプターデモ “scene-Helicopter-27.Blend” を用います。

Blenderデモ

以下の画像の場所からダウンロードできます。感謝!。

Blenderの設定

設定をエクスポートする前に、Blenderの設定を確認します。

1. RenderEngineの設定

使用するレンダリングエンジンを選択します。 ここは Cycles Render を選択します (ヘリコプターデモの場合は既にCycles Renderが選択されています)。

2. Dimensionsの設定

出力する画像の解像度を設定します。 Resolution の項目を100%にし、XYに画像の解像度を設定します。

3. Samplingの設定

Pattern で使用する乱数のサンプラーを選択し、 SettingsS に乱数のシード値を指定しします。 SamplesRender に1ピクセル当りのレイのサンプル数を指定します。

4. LightPathsの設定

様々な経路長のパスを考慮するため、Full Global Illumination を選択します。

5. PostProcessingの設定

画像処理前のレンダリング画像を取得するため、 Post Processingのチェックを全て外しておきます。

6. Surfaceの設定

Cyclesエンジンの初期設定で、うっすらとした環境光が設定されている場合があります。 意図しない光源が設定されていると比較の際に困るので、 BackgroundStrength0 になっているか確認しておきます。

7. Performanceの設定

パフォーマンスの比較を行う場合は、Threads で使用するスレッド数を指定します。

リファレンスの作成

以上の設定を確認してレンダリングを行います。

自作レンダラーとの比較のために、マテリアルの設定をもう少しシンプルにしておきます。

これで、設定をエクスポートする準備が整いました。

Blenderエクスポーターの作成

Blenderから設定をエクスポートするツールを作成します。

BlenderはPythonインタプリタを内蔵しており、 Blender内でPythonスクリプトを実行することができます。 そこで、Blenderから設定をエクスポートするPythonスクリプトを作成してみます。

Pythonスクリプトの実行

Pythonスクリプトは、Blender起動時に引数として渡して実行することができます。 試しにスクリプトを実行してみます。 以下は.Blendファイルのファイル名を表示するスクリプトです。

端末を開いて、Blenderに引数として渡して実行します。

% blender --background Helicopter.blend --python exporter1.py
Read blend: Helicopter.blend

## FileName: Helicopter

Blender quit

Blenderの実行ファイルの場所は、Windowsなら C:\Program Files\Blender Foundation\Blender\blender.exe、 macOSなら/Applications/Blender/blender.app/Contents/MacOS/blenderにあると思います。

--background オプションを指定すると、BlenderのGUIを起動せずに Pythonスクリプトを実行することができます。

実行する際のポイントは、Helicopter.blendよりも後に--python exporter1.py を指定することです。 Helicopter.blendの設定が読み込まれた後にスクリプトが実行されます。

スクリプト内でimport bpy を実行することで、BlenderのAPIを利用できるようになります。 利用できるAPIの詳細は Blender Documentation Python API で参照できます。

シーンデータの取得

.Blendファイルを読み込んだ際の各種設定は、 bpy.data に記録されます (bpy.dataの構造については BlendData を参照して下さい)。

各シーンの情報については bpy.data.scenes に記録されています。 (bpy.data.scenesの構造については Scene を参照して下さい)。

特定のシーンデータを取得する場合はシーン名を指定することで取得できます。 シーン名は Outliner で確認できます。

以下はSceneという名前のシーンから情報を取得するスクリプトの例です。

最初に、 scene = bpy.data.scenes[scene_name] で欲しいシーン情報を取得し、 シーンの各種情報を取得しています。 上のスクリプトでは、情報を標準出力に出力しているだけですが、 実際には自作レンダラーのフォーマットに加工してファイルに出力します。

実行結果は以下のようになります。

% blender --background Helicopter.blend --python exporter2.py
Read blend: Helicopter.blend

Image resolution: 1920 x 1080
Sampler: CORRELATED_MULTI_JITTER
Sampler seed: 123456789
Threads: 4
Render samples: 512

Blender quit

もし欲しい情報のAPIがわからない場合は、 BlenderのGUI上でその情報を編集するためのボタンやメニューの上に マウスカーソルを置いてみて下さい。 しばらくするとポップアップが現れて、その情報を利用するためのAPIを表示してくれます。

カメラ

シーンのアクティブカメラの情報は、scene.camera に入っています。

scene.camera 自体は Object 構造になっており、 カメラの位置や回転などの情報が入っています。 scene.camera.data にはカメラの画角やレンズなどの情報が入っています (Camera を参照して下さい)。

実行結果は以下のようになります。

% blender --background Helicopter.blend --python exporter3.py
Read blend: Helicopter.blend

X axis rotation: 4.169113636016846 radian.
Y axis rotation: 3.1490345001220703 radian.
Z axis rotation: 3.9584245681762695 radian.
position: (1.483912467956543, -1.3075470924377441, 1.6609470844268799)
Horizontal field of view: 0.8575560450553894
Vertical field of view: 0.5033799409866333

Blender quit

オブジェクト

シーン内のオブジェクトのデータは scene.objects に入っています (各オブジェクトの構造は Object を参照して下さい)。 例えば、シーン内のオブジェクト名を全て表示するスクリプトは、

のようになります。

試しに上画像のOutlinerのシーンで実行してみると、

blender --background TestScene.blend --python exporter4.py
Read blend: TestScene.blend

ObjectB
ObjectA
ObjectA-2
ObjectA-1
Lamp
Camera

Blender quit

scene.objects にはオブジェクト同士の親子関係に関わらず全てのオブジェクトが 入っています。 親子関係を考慮する場合は、

のように書くことができます。 obj.parent には親となるオブジェクトが入っており (トップレベルの場合は None)、 obj.children は子となるオブジェクトの配列(子がいない場合は要素を持たない配列) となっており、これらを用いて親子階層を意識した処理を書くことができます。

実行結果は

% blender --background TestScene.blend --python exporter5.py
Read blend: TestScene.blend
ObjectB
ObjectA
  child node: ObjectA-1
  child node: ObjectA-2
Lamp
Camera

Blender quit

となります。

次は、オブジェクト内の情報をエクスポートしてみます。

obj.type を見ることで、どんな種類のオブジェクトなのかがわかります (例えば ‘MESH’、‘CURVE’、’CAMERA’など。Type of Object)。 また、この種類の違いによって、 obj.data の構造も異なります (’MESH’なら Mesh 、’CURVE’なら Curve など)。

bpy.ops.export_scene.obj はオブジェクトを .obj 形式で書き出すAPIです (Export Scene Operators を参照)。 オブジェクトの種類毎に処理を変えるのは大変なので、 ここでは bpy.ops.export_scene.obj を用いて ’MESH’や’CURVE’などのオブジェクトは全て .obj 形式に書き出すようにしました。

以下は実行結果の一部を切り出したものです。

% blender --background Helicopter.blend --python exporter6.py
Read blend: Helicopter.blend

...

Bolt+Nut.011
  Material name: Metal
  Material index: 3
    (  0.0409 sec |   0.0000 sec) OBJ Export path: 'Bolt+Nut.011.obj'
          (  0.2900 sec |   0.2481 sec) Finished writing geometry of 'Bolt+Nut.011'.
      (  0.2964 sec |   0.2550 sec) Finished exporting geometry, now exporting materials
      (  0.2964 sec |   0.2550 sec) OBJ Export Finished

Progress: 100.00%

...

Blender quit

オブジェクト“Bolt+Nut.011”が“Bolt+Nut.011.obj”として書き出されました。 自作レンダラーでは、このobjファイルを読み込むように設定をエクスポートしています。

マテリアル

マテリアルの情報は bpy.data.materials に入っています (各マテリアルの構造は Material を参照して下さい)。 Blenderのマテリアルデータは複雑で、また、自作レンダラーのマテリアルデータとは 取るパラメータも異なるため、 マテリアルデータに関してはエクスポートは行わず後から手動で設定するようにしました。

Cyclesレンダラーと自作レンダラーの比較

エクスポートした設定を自作レンダラーに読み込ませて実際にレンダリングしてみました。

以下はCyclesレンダラーとの比較になります。

上画像がCyclesレンダラーでレンダリングした画像で、 下画像が自作レンダラーでレンダリングした画像です。 カメラ設定やオブジェクトの場所は一致させられたと思います。

しかし、レンダラー間でのマテリアルが取るパラメータの違いや係数の違い、 色空間設定の違いから、 結果を完全に一致させるにはまだ工夫が必要だと感じました。 一応今のままでも簡単な比較には使えそうです。

私のレンダラー向けのエクスポーターは、 blender_scene_exporter.py にあります。 参考にどうぞ。

参考

  1. Blender Documentation Python API: https://docs.blender.org/api/current/

  1. この記事内で使用したBlenderのバージョンは 2.79b です。

2017年12月17日日曜日

レンダリング方程式の変数変換

rendering_equation_jacobian

レンダリング方程式の変数変換

このページは レイトレ Advent Calendar 2017 の記事である. 本ページでは,レンダリング方程式の変数変換について考える.

記号について

本ページでは以下の記号は定義済として扱う.

記号 説明
\(\omega_{n}\) 単位法線ベクトル
\(\omega_{i}\) 入射方向の単位ベクトル
\(\omega_{o}\) 反射方向の単位ベクトル
\(L_{o}\) ある点からの反射光の放射輝度
\(L_{e}\) ある点が自発光している場合の放射輝度
\(L_{i}\) ある点への入射光の放射輝度
\(S^{2}\) 半球上の微小立体角の集合
\(M\) 空間中にある物体上の微小面積の集合
\(f_{r}\) BRDF

微小立体角を積分する形のレンダリング方程式は式1で表される.

\[ \begin{align} L_{o} \left( x, \boldsymbol{\omega_{o}} \right) = L_{e} \left( x, \boldsymbol{\omega_{o}} \right) + \int_{S^{2}} f_{r} \left( x, \boldsymbol{\omega_{i}}, \boldsymbol{\omega_{o}} \right) L_{i} \left( x, \boldsymbol{\omega_{i}} \right) cos \theta_{\omega_{n}} d \sigma \left( \boldsymbol{\omega_{i}} \right) \tag{1} \end{align} \]

式1では, 方向 \(\omega_{i}\) を中心とした微小立体角を半球上で積分することによって 放射輝度の計算を行っている.

一方で,入射方向 \(\omega_{i}\) ,出射方向 \(\omega_{o}\) の代わりに, 下図のように空間中の3点を接続する形のレンダリング方程式を考えることもできる.

3点接続のレンダリング方程式では, 点 \(x^{''}\) を中心とした微小面積を 空間中の物体上で積分することで放射輝度の計算を行う.

微小面積の積分では,点 \(x\) から半球を通して見える範囲の微小面積のみを積分したい. そこで次の関数を導入する.

\[ \begin{align} V \left( x \leftrightarrow x^{''} \right) = \begin{cases} 1 \ \ \left( x と x^{''} の間に遮蔽物が無い場合 \right) \\ 0 \ \ \left( x と x^{''} の間に遮蔽物がある場合 \right) \end{cases} \tag{2} \end{align} \]

この関数によって,点 \(x\) から半球を通して見えない微小面積をカットする.

次に,微小立体角 \(d \sigma\) と微小面積 \(dA\) の関係は,下図から式3のように表すことができる.

\[ \begin{align} d \sigma = \frac{dA cos \theta_{\omega_{n^{''}}}}{\left| x - x^{''} \right|^{2}} \tag{3} \end{align} \]

式2,3から,式1を基にして微小面積を積分する形のレンダリング方程式は

\[ \begin{align} L_{o} \left( x \rightarrow x^{'} \right) = L_{e} \left( x \rightarrow x^{'} \right) + \int_{M} f_{r} \left( x^{'} \rightarrow x \rightarrow x^{''} \right) L \left( x \leftarrow x^{''} \right) V \left( x \leftrightarrow x^{''} \right) \frac{cos \theta_{\omega_{n}} cos \theta_{\omega_{n^{''}}}} {\left| x - x^{''} \right|^{2}} dA \left( x^{''} \right) \tag{4} \end{align} \]

となる.

まとめ

本ページではレンダリング方程式の変数変換について簡単にまとめた. 積分の変数変換ができれば,論文の式の検証ができたり自分で新しい式を考えることができて 便利です.

参考文献

  1. James T. Kajiya: The rendering equation (1986)

2017年8月13日日曜日

Layered Diffuse Surface (Interfaced Lambertian BRDF)

Layered Diffuse Surface (Interfaced Lambertian BRDF)

このページは レイトレ合宿5‽ アドベントカレンダー第9週目の記事です.

本ページでは,モンテカルロレイトレーシング法において, 上図のような鏡面反射する層を持った拡散反射基板の表面 (Interfaced Lambertian BRDF) をレンダリングする方法を説明します. Interfaced Lambertian BRDF の実装に必要な パラメータや数式,重点的サンプリングの方法を説明します.

記号について

本ページでは以下の記号は定義済として扱います.

記号 説明
\(\omega_{n}\) Macrosurfaceの単位法線ベクトル
\(\omega_{m}\) Microsurfaceの単位法線ベクトル,本ページではマイクロファセット法線と呼ぶ
\(\omega_{i}\) 入射方向の単位ベクトル
\(\omega_{o}\) 反射方向の単位ベクトル
\(f_{r}\) BRDF
\(n_{i}\) 入射側の媒体の屈折率
\(n_{o}\) 透過側の媒体の屈折率
\(n\) \(\frac{n_{o}}{n_{i}}\)
\(F\) フレネル項

本モデルについて

Interfaced Lambertian BRDF は,下図のようにフレネル反射層と拡散反射基板から構成されます. このモデルは,フレネル反射層の境界での鏡面反射と 層内部の多重反射後の散乱を考慮しています. フレネル反射層と拡散反射基板はそれぞれ鏡面反射,拡散反射のMicrofacetモデルで表現され, 必要なパラメータもそれぞれのMicrofacetモデルのパラメータを組み合わせるだけの シンプルなモデルとなっています.

Interfaced Lambertian BRDF

Interfaced Lambertian BRDF の評価は, 入射光に対する鏡面反射 \(f_{s}\) と, 多重反射後の散乱 \(f_{b}\) に分けて行います.

\[ \begin{align} f_{r} \left( x, \omega_{i}, \omega_{o} \right) = f_{s} \left( x, \omega_{i}, \omega_{o} \right) + f_{b} \left( x, \omega_{i}, \omega_{o} \right) \tag{1} \end{align} \]

入射光に対する鏡面反射 \(f_{s}\) には, 鏡面のMicrofacetモデルを用います (Beckman,GGXなど) .

\[ \begin{align} f_{s} \left( x, \omega_{i}, \omega_{o} \right) = \frac{1}{4 \left( \omega_{i} \cdot \omega_{n} \right) \left( \omega_{o} \cdot \omega_{n} \right)} F \left( \omega_{i}, \omega_{m} \right) G_{2} \left( \omega_{i}, \omega_{o}, \omega_{m} \right) D \left( \omega_{m} \right) \tag{2} \end{align} \]

多重反射後の散乱 \(f_{b}\) は, 拡散反射面のMicrofacetモデルに, 入射光の内部への屈折 \(T \left( \omega_{i}, \omega_{m} \right)\) と 内部での多重反射 \(\frac{1}{\pi n^{2} \left( 1 - r_{i} \rho \right) }\) と, 内部から外部へ散乱する時の屈折 \(T \left( \omega_{o}, \omega_{m} \right)\) を 考慮したものになります.

\[ \begin{align} f_{b} \left( x, \omega_{i}, \omega_{o} \right) = & \frac{\rho}{\pi n^{2} \left( 1 - r_{i} \rho \right)} \int T \left( \omega_{i}, \omega_{m_{b}} \right) T \left( \omega_{i}, \omega_{m_{b}} \right) D \left( \omega_{m_{b}} \right) G \left( \omega_{i}, \omega_{o}, \omega_{m_{b}} \right) \frac{ \left| \omega_{i} \cdot \omega_{m_{b}} \right| \left| \omega_{o} \cdot \omega_{m_{b}} \right| } { \left| \omega_{i} \cdot \omega_{n} \right| \left| \omega_{o} \cdot \omega_{n} \right| } d \omega_{m_{b}} \\ r_{i} = & 1 - \frac{1 - r_{e}}{n^{2}} \\ r_{e} = & \frac{1}{2} - \frac{2 n^{3} \left( n^{2} + 2 n - 1 \right)} {\left( n^{2} + 1 \right) \left( n^{4} - 1 \right)} + \frac{\left( n - 1 \right) \left( 3 n + 1 \right)} {6 \left( n + 1 \right)^{2}} + \frac{8 n^{4} \left( n^{4} + 1 \right)} {\left( n^{2} + 1 \right) \left( n^{4} - 1 \right)^{2}} ln \left( n \right) + \frac{n^{2} \left( n^{2} - 1 \right)^{2}} {\left( n^{2} + 1 \right)^{3}} ln \left( \frac{n - 1}{n + 1} \right) \tag{3} \end{align} \]

この式には積分が含まれているため,モンテカルロ法を用いて Microfacet法線 \(\omega_{m_{b}}\) をサンプリングして解きます.

\[ \begin{align} f_{b} \left( x, \omega_{i}, \omega_{o} \right) = \frac{\rho}{\pi n^{2} \left( 1 - r_{i} \rho \right)} \frac{ T \left( \omega_{i}, \omega_{m_{b}} \right) T \left( \omega_{i}, \omega_{m_{b}} \right) D \left( \omega_{m_{b}} \right) G \left( \omega_{i}, \omega_{o}, \omega_{m_{b}} \right) \frac{ \left| \omega_{i} \cdot \omega_{m_{b}} \right| \left| \omega_{o} \cdot \omega_{m_{b}} \right| } { \left| \omega_{i} \cdot \omega_{n} \right| \left| \omega_{o} \cdot \omega_{n} \right| }} {p \left( \omega_{m_{b}} \right) } \tag{4} \end{align} \]

重点的サンプリング (Importance Sampling)

Interfaced Lambertian BRDF の重点的サンプリングは, フレネル反射層の反射分布,もしくは拡散反射基板の反射分布に従ってサンプリングを行います.

まず最初は,フレネル反射層,拡散反射基板どちらの反射分布に従ってサンプリングを行うかを決定します. ここでは,フレネル反射層,拡散反射基板それぞれの総反射率をもとに確率的に決定します.

\[ \begin{align} R_{s} = & r_{e} \\ R_{b} = & \frac{ \left( 1 - r_{e} \right) ^{2} \rho } { n^{2} \left( 1 - r_{i} \rho \right) } \tag{5} \end{align} \]

\(R_{s}\)\(R_{b}\) はそれぞれフレネル反射層,拡散反射基板の総反射率です. \(u_{1} \in [0, 1)\) となる一様乱数を用いて, \(u_{1} < \frac{R_{s}}{ \left( R_{s} + R_{b} \right) }\) となる場合は フレネル反射層の反射分布に従ったサンプリングを行い,そうでない場合は, 拡散反射基板の反射分布に従ったサンプリングを行います.

フレネル反射層の反射分布のpdf \(p_{s} \left( \omega_{o} \right)\) や 反射方向のサンプリング方法は層のMicrofacetモデルに従う (Beckman,GGXなど). 反射方向のサンプリングの方法については,アドベントカレンダー第8回Rough Specular Surface (GGX BRDF) を参照してください.

拡散反射基板の反射分布のpdf \(p_{b} \left( \omega_{o} \right)\) は,

\[ \begin{align} p_{b} \left( \omega_{o} \right) = \frac{ \left( \omega_{o} \cdot \omega_{n} \right) }{\pi} \tag{6} \end{align} \]

であり, \(p_{b} \left( \omega_{o} \right)\) に従って反射方向をサンプリングします. 反射方向のサンプリング方法は Smooth Diffuse Surface (Lambert BRDF) を参照してください.

最終的な反射方向のpdf \(p \left( \omega_{o} \right)\) は,

\[ \begin{align} p \left( \omega_{o} \right) = \frac{R_{s}}{ \left( R_{s} + R_{b} \right) } p_{s} \left( \omega_{o} \right) + \frac{R_{b}}{ \left( R_{s} + R_{b} \right) } p_{b} \left( \omega_{o} \right) \tag{7} \end{align} \]

となります.

まとめ

パラメータ

名称 記号
拡散反射基板の反射率 \(\rho\) \(0 \le \rho \le 1\)
表面粗さ \(roughness\) \(0 < roughness \le 1\)

レンダリング

名称 記号
BRDF \(f_{r} \left( x, \omega_{i}, \omega_{o} \right)\) \(f_{s} \left( x, \omega_{i}, \omega_{o} \right) + f_{b} \left( x, \omega_{i}, \omega_{o} \right)\)
pdf \(p \left( \omega_{o} \right)\) \(\frac{R_{s}}{ \left( R_{s} + R_{b} \right) } p_{s} \left( \omega_{o} \right) + \frac{R_{b}}{ \left( R_{s} + R_{b} \right) } p_{b} \left( \omega_{o} \right)\)
反射レイのWeight \(\frac{f_{r} \left( x, \omega_{i}, \omega_{o} \right) \left( \omega_{o} \cdot \omega_{n} \right)} {p \left( \omega_{o} \right)}\) \(\frac{f_{r} \left( x, \omega_{i}, \omega_{o} \right) \left( \omega_{o} \cdot \omega_{n} \right)} {p \left( \omega_{o} \right)}\) (特にキャンセルできる項は無い)

以上,Interfaced Lambertian BRDFについて簡単にまとめました.

参考文献

  1. Daniel Meneveaux: Rendering Rough Opaque Materials with Interfaced Lambertian Microfacets (2017)

2016年8月7日日曜日

BRDFのテスト

BRDFのテスト

このページは レイトレ合宿4!? アドベントカレンダー第7週目の記事です. レイトレ合宿まで1ヶ月をきりました. そろそろ提出用のシーン作成やレンダラーのデバッグに集中しているころではないでしょうか.

私はレンダラーの実装に不具合があるかどうかは, レンダリングした画像を見て違和感を持つかどうかで判断しています. この方法は主観的であり, 画像に違和感を持たなくても不具合が含まれている可能性は大いにあります. そこで,自動テストを実装して客観的に不具合の判定をするようにします.

このページでは物理ベースなBRDFの自動テストについて考えます.

物理ベースなBRDF

初めに,物理ベースなBRDFについて確認しておきます.

定義済み記号 説明
\(\omega_{n}\) 法線ベクトル
\(\omega_{i}\) 入射方向の単位ベクトル
\(\omega_{o}\) 反射方向の単位ベクトル
\(f_{r}\) BRDF

物理ベースなBRDFは,反射する放射エネルギーが負になることはないので

\[ \begin{align} \forall \boldsymbol{\omega_{i}} \ \forall \boldsymbol{\omega_{o}}, \ f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ cos \theta_{\boldsymbol{\omega_{o}}} \ge 0 \end{align} \]

\(cos \theta_{\boldsymbol{\omega_{o}}} \ge 0\) であるため,

\[ \begin{align} \forall \boldsymbol{\omega_{i}} \ \forall \boldsymbol{\omega_{o}}, \ f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ge 0 \tag{1} \end{align} \]

となります. また,物理ベースなBRDFは以下の2つの特性を持ちます.

ヘルムホルツの相反性 (Helmholtz Reciprocity)

BRDFの値は入射・出射方向が入れ替わっても変化しません.

\[ \begin{align} f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) = f_{r}(x, \boldsymbol{\omega_{i}} \leftarrow \boldsymbol{\omega_{o}}) \tag{2} \end{align} \]

エネルギー保存則 (Energy Conservation)

ある点から反射した総出射光束は総入射光束以下となります.

\[ \begin{align} \forall \boldsymbol{\omega_{i}}, \ \int_{\Omega} f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ cos \theta_{\boldsymbol{\omega_{o}}} d \sigma(\boldsymbol{\omega_{o}}) \le 1 \tag{3} \end{align} \]

モンテカルロ法によるテスト

\((1)\)\((2)\)\((3)\)の自動テストを実装します. 式\((3)\)は解析的に解くことが難しいため,モンテカルロ法を用いてテストを行います. モンテカルロ法を用いた場合,式\((3)\)は以下のように書きかえられます.

\[ \begin{align} \forall \boldsymbol{\omega_{i}}, \frac{1}{N} \sum_{k=1}^{N} \frac{f_{r}(x, \boldsymbol{\omega_{ok}} \leftarrow \boldsymbol{\omega_{i}}) cos \theta_{\boldsymbol{\omega_{ok}}}} {p_{u} \left( \boldsymbol{\omega_{ok}} \right)} \le 1 \tag{4} \end{align} \]

ここで,\({p_{u} \left( \boldsymbol{\omega_{o}} \right)}\)\(\boldsymbol{\omega_{o}}\)の確率密度関数 (PDF) です. このテストでは,半球上の方向を一様にサンプルするようにします \(\left({p_{u} \left( \boldsymbol{\omega_{o}} \right)} = \frac{1}{2 \pi} \right)\)

C++言語で GoogleTest を用いたテストの擬似コードを記述します.

// エネルギー保存則のテスト
TEST(BrdfTest, EnergyConservationTest)
{
  Brdf brdf; // テスト対象であるBRDFの宣言,初期化
  const int N = 1000000;
  // 様々な入射方向に対してテストを行う
  for (int i = 0; i < N; ++i) {
    const float energySum = 0.0;
    const Vector3 vin = sampleDirectionUniformly(); // 半球上の方向を一様にサンプル
    for (int o = 0; o < N; ++i) {
      const Vector3 vout = sampleDirectionUniformly();
      const float f = brdf.f(vin, vout); // BRDFの評価
      ASSERT_LE(0.0, f); // 式(1)のテスト
      const float cosTheta = dot(normal, vout);
      energySum += f * cosTheta * (2.0 * pi); // pi は円周率
    }
    const float energy = energySum / (float)N;
    ASSERT_GE(1.0, energy); // 式(4)のテスト
  }
}
// ヘルムホルツの相反性のテスト
TEST(BrdfTest, HelmHoltzReciprocityTest)
{
  Brdf brdf;
  const int N = 1000000;
  for (int i = 0; i < N; ++i) {
    const Vector3 vin = sampleDirectionUniformly();
    const Vector3 vout = sampleDirectionUniformly();
    const float f1 = brdf.f(vin, vout); // BRDFの評価
    const float f2 = brdf.f(-vout, -vin); // 方向を入れ替えたBRDFの評価
    const float error = 1.0e-5; // 許容できる誤差を定義
    ASSERT_FLOAT_EQ(f1, f2, error); // 式(2)のテスト
  }
}

実際にテストを実装する場合は, モンテカルロ法による\(N\)回試行の和を求める際に情報落ちの危険があるため, 補正加算も考慮したほうがいいと思います.

重点的サンプリング (Importance sampling) のテスト

BRDFによっては,反射方向を重点的サンプリングする式が定義してあることがあります. 重点的サンプリングの際の反射方向のPDFを \(p \left( \boldsymbol{\omega_{o}} \right)\) とします. 確率密度関数の性質から,

\[ \begin{align} \int_{\Omega} p \left( \boldsymbol{\omega_{o}} \right) d \sigma \left( \boldsymbol{\omega_{o}} \right) = 1 \tag{5} \end{align} \]

となります.式\((3)\)と同様に,解析的に解くことが難しいためモンテカルロ法を 用いて以下のように書きかえます.

\[ \begin{align} \frac{1}{N} \sum_{k=1}^{N} \frac{p \left( \boldsymbol{\omega_{ok}} \right)} {p_{u} \left( \boldsymbol{\omega_{ok}} \right)} \approx 1 \tag{6} \end{align} \]

擬似コードを以下に記述します.

// PDFのテスト
TEST(BrdfTest, PdfTest)
{
  Brdf brdf;
  const int N = 1000000;
  // 様々な入射方向に対してテストを行う
  for (int i = 0; i < N; ++i) {
    float pdfSum = 0.0;
    const Vector3 vin = sampleDirectionUniformly();
    for (int o = 0; o < N; ++o) {
      const Vector3 vout = sampleDirectionUniformly();
      const float pdf = brdf.pdf(vin, vout); // 反射方向のPDFを評価
      pdfSum += pdf * (2.0 * pi);
    }
    const float pdf = pdfSum / (float)N;
    const float error = 1.0e-5; // 許容できる誤差を定義
    ASSERT_FLOAT_NEAR(1.0, pdf, error); // 式(6)のテスト
  }
}

終わりに

BRDFのテストについて考えました. このテストで確認できたことは,物理に基づいたBRDFの定義を満たしているか ということであり,BRDFの反射分布が正しいかというところまでは確認できていません. ここからテスト項目を増やしていって,確認できることを広げていこうと思います.