C++ でレイマーチングをする

はじめに

この記事は C++ Advent Calendar 2019 の 1 日目の記事です。

最近、 Twitter で "C++ でレイマーチングをする" という話題が時々まわってきます。実際に自分がやるならこうやるかなあ的な内容で記事を書いてみようかと思います。

実際のところレイマーチング (プロシージャルグラフィックス) 自体より、それを実現するための

• 数学ライブラリ
• 並列実行
• 画面表示

についての紹介、解説記事になります。

サンプルプロジェクトはこちらです。

• aosoft/CppRayMarching

Windows 10 + Visual Studio 2019 ですが、実装したコードに関しては Windows に依存しないように書いたはずなので、他環境でも (依存ライブラリが対応していたら) 動作するはずです。

ビルド済バイナリ (Windows x64) も作っておきましたので、よろしければお試しください。コマンドラインオプションで数字 (0～3) を指定するとそれに応じた関数でレンダリングをします。指定なしは 0 です。 ESC キーで中断できます。 CPU 使用率は非常に高く、大抵の PC でほぼ 100% に達すると思われますのでご注意を。

• 2019/11/28 版

グラフィックス生成関数を実装する

OpenGL Mathematics (数学ライブラリ)

コードで映像を生成する場合、その実装は数学処理の塊なので、型 (ベクトル、行列) 定義や数学関数各種の実装が必要になります。これを独自に実装するから始めるのもよいですが、私は OpenGL Mathematics (GLM) を利用するのを強くお勧めします。

• OpenGL Mathematics
• g-truc/glm

GLM は GLSL との互換性を強く意識した設計で

• GLSL と同名のベクトル、行列型
• GLSL と同じ数学関数
• 各種便利数学関数 (perspective など)

となっていて OpenGL と合わせて使うととても親和性が高いですし、 GLSL で実装されたシェーダーを比較的容易に C++ へ移植できます。

GLM は GLSL との互換性が高いだけで、 OpenGL とは全く無関係のライブラリです。なので例えば Direct3D と合わせて使うことも問題なく可能です。実際、私は Direct3D と OpenGL とのクロス開発では GLM を利用して主要処理の共有化をしています。

Direct3D 向けには Windows SDK に含まれている "DirectX Math" があります。 DirectX Math も DirectX から独立した数学ライブラリですのでこちらを利用しても実現可能とは思いますが、実質 Windows 依存になってしまうのと、 GLM を使う場合と比較して自然な書き方にならないように思うのであまりお勧めはできません。

実際に実装する

GLM を include する際、その前に #define GLM_FORCE_SWIZZLE を定義します。これを定義しないと Swizzle が使えません。

#define GLM_FORCE_SWIZZLE

#include <glm/glm.hpp>

using namespace glm;

GLM は全て glm 名前空間下に定義されていますが、 GLSL からの移植の場合はいちいちつけるのが面倒なので using namespace をしましょう (.cpp で局所的に使うように) 。

レイマーチング関数は後述するレイマーチング実行処理部に関数ポインターで渡すことにします。今回のサンプルでは定義は次のようにしています。

using FnKernel = void (*)(glm::vec4& gl_FragColor, const glm::vec4& gl_FragCoord, float time, const glm::vec2 resolution);

Fragment Shader で必須な変数、定数や GLSL でレイマーチングをする際のお約束的定数を引数で渡すことで、内部実装はほぼ共通で記述可能になります。

GLSL からの移植の場合は次のような感じにするのが手っ取り早いです。

1. とりあえずコピペする
2. main 関数をエントリー定義に合わせる
3. コンパイルしてみてエラーを一つずつつぶす
   • 浮動小数点型の定数は C++ の場合 double になってしまうので接尾辞をつけて float にする (1.0 → 1.0f)
   • Swizzle は operator でのアクセスにしないといけない場合があるので、必要に応じて対応する (p.xy → p.xy())

GLM は原則 float での演算になりますが、演算子オーバーロードなどで float じゃないとうまくいかない (double から暗黙のキャストがきかない) 事があるので適宜 float にしましょう。

Swizzle は swizzle 型で公開されますが、原則として swizzle 型では入力型として扱えません。 vec 型に変換する operator () オーバーロードが定義されているので、 Swizzle を入力ソースとして使う場合は後ろに () を付けると解消します。

例として Sphere のレイマーチングの C++ (左) と GLSL (右) の比較です。簡単なコードではありますがほとんど差分がない事がわかるかと思います。

レンダリングを実行する

Fragment Shader のような処理はピクセル単位での独立した処理なので完全に並列実行ができます。よって CPU で実行する場合でもマルチコアを生かして並列実行をした方がよいです。

並列実行を実装するにはコンパイラの機能 (OpenMP) やライブラリ (Microsoft の PPL, Intel の TBB など) を利用するなどが考えられます。

シングルスレッド (基本)

いきなりスレッド分割することは普通はしないので、とりあえず単純にループをまわして実行します。

void RunKernelSingle(std::int32_t width, std::int32_t height, void* p, std::int32_t stride, float time, FnKernel fnKernel)
{
    auto resolution = glm::vec2(static_cast<float>(width), static_cast<float>(height));

    for (std::int32_t y = 0; y < height; y++)
    {
        auto p2 = reinterpret_cast<glm::vec4*>(static_cast<std::uint8_t*>(p) + (height - y - 1) * stride);
        for (std::int32_t x = 0; x < width; x++)
        {
            fnKernel(*p2, glm::vec4(static_cast<float>(x), static_cast<float>(y), 0.0f, 0.0f), time, resolution);
            p2++;
        }
    }
}

OpenMP

OpenMP はコンパイラの機能で対応コンパイラのコンパイルオプションで OpenMP を指定すると有効化されます。並列化したいループの前にちょこっと構文を付けると自動でいい感じにやってくれます。

#pragma omp parallel for
for (std::int32_t y = 0; y < height; y++)

OpenMP が有効になっていなくてもコンパイル自体は成功しますが生成されたコードは並列化されずそのままシングルスレッド動作になるので、その点は注意してください。

std::async で自前スレッド分割をする

自前で並列実行を実装する場合、大ざっぱには

• 必要数のスレッドを生成
• 生成したスレッドに応じて適切にタスク分割
• 分割したタスクをそれぞれのスレッドに渡して実行

のような感じになります。

ここで問題になってくるのが "スレッド生成" なのですが、一般的にスレッド生成はコストが高いため、処理毎に都度スレッド生成をしていると効率が悪くなります。なのでスレッドプールを用意してスレッドを再利用するのが常套的な手段になりますが、 C++ では標準ライブラリでスレッドプールがないようで、ここでまた困ってしまいました。自前で実装してもいいのですが微妙に面倒ですし流石に本題からズレすぎかなあと (boost 使えば、というのは見なかった事に) 。

いろいろ調べたり試した結果、 std::async が実装系依存ですがスレッドプールを使った実装が期待できることと、 Visual C++ では PPL で実際にスレッドプールで動作していたので、今回は std:async で並列実行のコードを書くことにしました。

std::atomic_int32_t y = 0;
for (size_t i = 0; i < parallelCount; i++)
{
    futures[i] = std::async(std::launch::async, [&y, width, height, p, stride, time, fnKernel]()
        {
            auto resolution = glm::vec2(width, height);
            while (true)
            {
                auto yy = y++;
                if (yy >= height)
                {
                    break;
                }
                auto p2 = reinterpret_cast<glm::vec4*>(static_cast<std::uint8_t*>(p) + (height - yy - 1) * stride);
                for (std::int32_t x = 0; x < width; x++)
                {
                    fnKernel(*p2, glm::vec4(static_cast<float>(x), static_cast<float>(yy), 0.0f, 0.0f), time, resolution);
                    p2++;
                }
            }
        });
}

分割はあまり粒度を細かくしすぎると同期やメモリアクセスの乱れなどでパフォーマンスに影響が出る可能性がありますが、ここではライン単位で分割するようにしました。 y を atomic 変数とし、 Y 座標を重複しないように取り合うようにし、それぞれのスレッドでライン別に処理するようにしています。横方向のメモリは連続していますが縦方向は連続していない (するかもしれないがする前提にはできない) ので。 y が atomic 変数で若干コスト的に気になるところではありますが、計測した感じだと特に支障もなく、コードもシンプルに書けているので今回はこのようにしました。

同時実行数 (parallelCount) は std::thread::hardware_concurrency() を参照しつつも最大数は適当な数で決め打ちにしています。

OpenCV の並列実行機能を利用する

PPL は VC++ 依存ですし、 TBB は別途ライブラリの導入が必要。今回は後段で OpenCV を使うので、 OpenCV の parallel_for での並列実行もしてみました。

cv::parallel_for_(cv::Range(0, height), [width, height, p, stride, time, fnKernel](const cv::Range& range)
    {
        auto resolution = glm::vec2(width, height);
        for (std::int32_t y = range.start; y < range.end; y++)
        {
            auto p2 = reinterpret_cast<glm::vec4*>(static_cast<std::uint8_t*>(p) + (height - y - 1) * stride);
            for (std::int32_t x = 0; x < width; x++)
            {
                fnKernel(*p2, glm::vec4(static_cast<float>(x), static_cast<float>(y), 0.0f, 0.0f), time, resolution);
                p2++;
            }
        }
    });

一般的な parallel for な定義になっていて、指定した並列実行したい範囲 (ここでは Y 方向の範囲) を OpenCV 側で適当に分割して分割単位の Range を渡してくれるので、その範囲毎に並列実行します。こちらの場合、分割ブロックが完全に独立しているので完了するまで衝突する要素がない反面、ブロック毎の処理の負荷が偏ると先に終了したタスクの CPU コアが遊んでしまう事が考えられます (大きな差異が出るほど偏りが出ることはまれだとは思いますが) 。

Windows における制限

Windows では "プロセッサーグループ" という仕組みにより 64 論理コアを超えるシステムを利用すると 1 プロセスに割り当てられる最大コア数が制限されるようになります。

• 64 コアを超える Windows 環境でマルチスレッド・プログラミングをしてみる

この制限を超えてフルに CPU をまわせるようにするには複数プロセスを生成してプロセス間をまたいでタスクをまわすようにしないといけないようです。

画面表示

レンダリングした結果をどう確認するかですが、ファイルに保存して別途ビューアーで見る、は正直面倒くさいので直接画面に表示したいところです。

画面表示となると OS の API か何かしらの GUI フレームワークを使う、とかなってしまいますが、画像処理をしたビットマップを単純に表示したい、という目的だと OpenCV が手軽でお勧めです。 OpenCV には画像処理の確認のために簡単な GUI 機能が実装されています。今回は OpenCV の本体機能である画像処理は一切使わずに便利機能である GUI 機能だけ利用します。

int main(void)
{
    cv::Mat img(768, 1024, CV_32FC4);

    cv::namedWindow("Window", cv::WINDOW_AUTOSIZE | cv::WINDOW_FREERATIO);

    Timer timer;
    while (cv::waitKey(1) < 0)
    {
        RunKernelParallel(img.cols, img.rows, img.ptr(), img.step, timer.Current(), Sphere);
        cv::imshow("Window", img);
    }
}

• namedWindow 関数でウィンドウ表示をする
• Mat クラスでイメージバッファを用意し、そのバッファに対してレンダリング関数で描画処理を実行する
• imshow 関数でウィンドウにイメージを表示する
• waitkey 関数でキー入力があるまで繰り返す

レイマーチング関数に限らず、とりあえず画像処理してその結果を見たい (描画パフォーマンスは無視) 場合はとても簡単便利でよいです。 OpenCV は非常に巨大なライブラリですが、使いたいものだけつまみ食い的に使えるようになっているのでこういった分りやすいところから入っていくのもよいかと思います。

ピクセルフォーマット

今回は GLSL のレイマーチングを C++ に移植する、というスタイルで進めていますが GLSL の場合は Fragment Shader が出力するピクセルはレイアウトが RGBA ですが OpenCV では BGRA 並びになっています。よって OpenCV の Mat に書き出す際に適切に変換する必要があります。

Mat はレンダリングの出力に合わせて 32bit float にしています。おそらくは 8bit uint で出力した方がよいのではと思いますが、画面表示は今回はそこまで速度追及する部分ではないので気にしないことにしました。サンプルプロジェクトでの時間計測はレンダリング処理のみとして画面表示部分は除外しています。

SIMD 最適化を有効化してみる

GLM を include する前に #define GLM_FORCE_INTRINSICS を定義しておくと GLM が SIMD 最適化版での動作になります。

• 2.11. GLM_FORCE_INTRINSICS: Using SIMD optimizations

#define GLM_FORCE_SWIZZLE
#define GLM_FORCE_INTRINSICS

#include <glm/glm.hpp>

通常はコンパイラの最適化オプション (VC++ だったら /arch) に準じた CPU 命令が使われるようですが、 GLM のプリプロセッサ定義を変えることでコンパイラ定義と関係なく、他の CPU 命令を利用することもできるようです。

例えば AVX2 を有効にしてコンパイルした場合、そのバイナリは AVX2 が使える PC じゃないと実行時にエラーで落ちるので、一般に頒布するバイナリで使う場合は配慮が必要です。また、必ずしも良い効果が得られるとは限らないので試して確認してください。

実装例

@kaneta1992 さんに許可をいただいて、 Shadertoy で公開されている GLSL のコードを C++ に移植してみました。

• Kaleidoscope Tunnel → C++
• Hologram Boxes → C++

単純な球体くらいだったらたいしたことないのですが、流石にこのクラスだと激烈に重いですね。

"Kaleidoscope Tunnel" はレイマーチングではないですが、入れさせてもらいました。

サンプルプロジェクトについて

ビルドの仕方

Windows 環境では下記手順でビルドができます (おそらく) 。

1. プロジェクトをチェックアウト
2. GLM の submodule を更新
3. プロジェクトのルートフォルダー直下に OpenCV のダウンロードサイト からダウンロードしたアーカイブの中身 (opencv フォルダー以下) をコピー、もしくはしジャンクション (シンボリックリンク) をはる。
   
4. CppRayMarching.sln を開いてビルド
5. ビルド出力先に OpenCV の DLL バイナリをコピー (Debug = opencv_world???d.dll / Release = opencv_world???.dll)

実行

コマンドラインオプションで数字 (0～3) を指定するとそれに応じた関数でレンダリングをします。映像サイズは 640x480 です。大きくすればリニアに重くなっていきます。

Sphere

Color Box

Kaleidoscope Tunnel

Hologram Boxes

.@kanetaaaaa さんのGLSLのコード(Hologram Boxes)をC++に移植してCPU(i7 7700HQ 4C/8T)で動かしてみると640x480で1枚2秒くらいかかりました。 pic.twitter.com/SCHlzEOqM0

— TAN-Y(たに) (@TANY_FMPMD) November 28, 2019

動作設定の変更

本当はコマンドラインで設定できるようにすべきなのですが、面倒だったのでソースコードを変更してビルドし直してください。

• defines.h でプリプロセッサによる挙動変更ができます
• レンダリングサイズは main 文の先頭の方に定数定義で設定できるようにしています

おわりに

CPU でレイマーチングは実際のところパフォーマンスが悪すぎて使いどころが限定されそうな気はします。

CPU は確かに GPU に比べて遅いですが、極端なペナルティやメモリ参照の制約もほとんどありませんし、 1 関数の実行時間に制限もありません (シェーダーは長すぎるとドライバーでエラーになる) 。今回は GLSL からの移植というスタンスでやっていますが、柔軟性という事では CPU の方が有利ですので、 CPU ならではの使い方を考えていけばよいのではないかと思いますし、今回の記事で挙げた実現するための各技術は汎用的に使えるものと思いますので知っておいて損はないと思います。

個人的に OpenGL Mathematics はものすごく便利だと思うので、 C++ でグラフィックス処理をやられる場合は是非使ってみてください。

今回の記事作成にあたり @kaneta1992 さんの GLSL コードを利用させていただきました。また、 @doxas さんのたくさんの記事も参考にさせていただきました。ありがとうございました。