ShirotsuのCG学習録

ShirotsuのCG学習録

オフライン/リアルタイムレンダリング, C++などの学習記録

CSゲーム業界で働きたい理系大学生へ捧げる本たち

これから何を学べばいいか困っている大学生に向け、独断と偏見でオススメ書籍をまとめた記事です。主に「授業でC言語とかやって、自分でC++もちょっとずつ触り始めたけど、どうすれば...」という方に読んで欲しいです。

ちょっとタイトル詐欺ですが、「自分が何を読んだか」という備忘録的な意味も込めています。

また、Shirotsu自身しっかり読んでいない本も少々混在しています。公開することに意義があると思うので、どうかご了承ください...

(2020/09/26 追記)
Twitterで予想以上の反響があり、「この本もどうか」というリプライを方々からいただきました。
せっかくなのでここに順次追加していきます。が、Shirotsu自身が完全に内容を知らないためにまだ内容を保証できないという意味で [未読]タグをつけておきます。

(2020/09/28 追記)
書籍の量が多くなってきたので、カテゴリの細分化とリネームを行ないました。

(2020/09/28 追記)
これ以降この記事の更新はしません。
代わりに、この記事の改良版を作りました。
そちらの方が項目の数が多くなっていますので、ぜひご覧ください。

入門レベル ～何を学べばいいか分からない段階～

中級レベル ～ある程度慣れてきた段階～

上級レベル ～だいぶ慣れて自信もついてきたら～

興味に応じて読むべき書籍

入門レベル

中級レベル

上級レベル

  • GPU Gems 日本語版 & GPU Gems 2 日本語版 & GPU Gems 3 日本語版 [未読]
    シェーダ含む、GPUでのあらゆるグラフィックス計算の術が記されたバイブル的名著。鈍器。
    1~3の続き物ですが、日本語版は1が絶版です。また販売価格も高く手に取りにくいです。大学図書館に置いてあったら手を出してみるくらいの気持ちでよいでしょう。
    ちなみに、こちらも1~3すべて英語版が無料公開されています。

辞書的に使える書籍

入門レベル

中級レベル

上級レベル

  • リアルタイムレンダリング 第4版
    前述のCG-ARTSのCG本の上級向けバージョンといった内容で、最近のトレンド的な話も多くあります。この本の内容をほぼ理解していれば神になれると僕は思います。ただ、翻訳の質があまりよくないのが玉に瑕。

Embree3 入門解説 処理フロー概要

※ この記事は、筆者が以前Qiitaに投稿したものの再掲版です。 ※

はじめに

Embreeとは、Intelが提供する高速なレイトレーシングを実現するC言語(C99)のライブラリです。
EmbreeはCPUのみで処理を実行し、非常に高いパフォーマンスを発揮してくれます。
この記事では、これからEmbree 3.x系を利用する人へ向け、C++による実装を処理の流れに沿って解説します。

また、この記事はレイトレ Advent Calendar 2019の記事です。
興味がある方はそちらもどうぞ！

想定読者

この記事の想定読者は以下の通りです。関連するトピックの用語を前置き無く使用しますのでご了承ください。

  • Embreeを利用し始めたいが勝手が分からず困っている方
  • 基本的なC++の記法が分かる方
  • レイトレーシングの基本的な描画フローが分かる方
  • Graphics PipelineのVertex Buffer, Index Bufferを知っている方

開発環境

この記事上のサンプルコードを動作させた環境です。

Embree APIバージョン

Embreeは、2018年3月に発表された2.x系から3.0へのアップデートの際にAPIが刷新されました。
2.x系と3.x系では処理の流れは概ね同じですが、2.x系のAPIのほとんどが3.x系で利用できません。
この記事では、Embree 3.x系に沿って解説します。(以降Embree3と表記)
筆者も大いに参考にさせていただいたEmbree 2.x系の解説記事はこちら。

描画モデル

tetra_bottom

今回はこのモデル(頂点数6, ポリゴン数8)を利用して解説します。
以下にC++コードで便宜的に示します。

vertex[0] = { -0.5f,  0.0f,  0.0f };
vertex[1] = { -0.0f,  1.0f,  0.0f };
vertex[2] = {  0.5f,  0.0f,  0.0f };
vertex[3] = {  0.0f, -1.0f,  0.0f };
vertex[4] = {  0.0f,  0.0f,  0.5f };
vertex[5] = {  0.0f,  0.0f, -0.5f };
index[0] = { 1, 0, 4 };
index[1] = { 1, 4, 2 };
index[2] = { 3, 0, 4 };
index[3] = { 3, 4, 2 };
index[4] = { 1, 5, 0 };
index[5] = { 1, 2, 5 };
index[6] = { 3, 5, 0 };
index[7] = { 3, 2, 5 };

環境構築

公式ページから依存ファイルをダウンロードし、お好きなディレクトリへインストールしてください。
ヘッダファイル, リンクファイルの設定もお忘れなく。

ちなみに、旧バージョンは公式開発リポジトリのReleasesから入手できます。

Embree3のオブジェクト

Embree3にはDevice, Scene, Geometry, Buffer, BVH というオブジェクトがあります。
(これらはEmbree3で特殊な意味を持つので、この記事上では曖昧さ回避のため便宜的に「Embree3オブジェクト」と表記します)
これらに対応する型のインスタンスは通常rtcNewXXX(XXXにDevice等が入る)という関数で生成し、その返り値のハンドルをユーザが利用します。
また、これらのオブジェクトは参照カウント方式で管理できます。rtcNewXXXで生成した直後はカウントが1で、rtcRetainXXX関数でインクリメント、rtcReleaseXX関数でデクリメントします。リリース時にカウントが0ならばインスタンスが解放されます。
(この一連の操作を「ユーザ側で管理」と表現します)

C++でラッパークラスを作るならばコンストラクタとデストラクタで操作すると使い勝手が良さそうです。

処理フロー概要

embree_qiita_api

Embree3ではこのような流れでプログラムを記述していきます。
下記の各解説では、利用する構造体と関数シグネチャ -> 呼び出し側のコード例 という順で表記します。

// 宣言されているヘッダファイル名
struct Hoge;

Hoge CreateHoge();

また、「Embree3のKernel」という単語を「APIコールによってEmbree3の内部で処理をする部分」という意味で用います。

デバイス生成

Embree3オブジェクトの一つ、デバイスを生成します。
返り値はハンドルとして受け取り、ユーザ側で管理します。

// embree3/rtcore_device.h
struct RTCDevice;
RTCDevice rtcNewDevice(const char* config);

RTCDevice device = rtcNewDevice(nullptr);

このデバイスは、他のEmbree3オブジェクトの生成に利用します。
また、引数に特定の文字列を渡すことでconfiguration)ができます。
NULLやnullptrを渡すとデフォルト設定になります。

RTCDevice device = rtcNewDevice("threads=1,isa=avx"); // コンマ区切りで複数指定

他にもデバイス生成時にconfiguration文字列を渡す方法があります。
数字の大きい方が優先され、小さい方は無視されます。

  1. rtcNewDeviceの引数に渡した文字列
  2. アプリケーションと同じディレクトリにある.embree3拡張子のファイル
  3. ホームディレクトリにある.embree3拡張子のファイル

具体的な設定項目については割愛させていただきます。
公式リファレンスのrtcNewDeviceの項に列挙されていますので、気になる方は覗いてみて下さい。
Embree3のKernelが利用するSIMD命令セットアーキテクチャやスレッド数の設定などができます。

シーン初期化

Embree3オブジェクトの一つ、シーンを生成します。
これは複数の「3D空間上のオブジェクト」を保持します。
返り値はハンドル(RTCScene型)として受け取り、ユーザ側で管理します。

// embree3/rtcore_scene.h
struct RTCScene;
RTCScene rtcNewScene(RTCDevice device);

RTCScene scene = rtcNewScene(device);

後に scene に「3D空間上のオブジェクト」を1つ以上アタッチし、レイトレースで利用します。

ジオメトリ初期化

ジオメトリ生成

「3D空間上のオブジェクト」ことジオメトリを生成します。これもEmbree3オブジェクトです。
返り値のハンドル(RTCGeometry型)を一時的に受け取りますが、この後の「シーンへのジオメトリ登録」の直後に解放し、ユーザ管理から切り離すことが できます。(理由は後述)

// embree3/rtcore_geometry.h
struct RTCGeometry;
RTCGeometry rtcNewGeometry(RTCDevice device, RTCGeometryType type);

RTCGeometry geometry = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_TRIANGLE);

rtcNewGeometryの第2引数にはEmbree3が定義しているenum値RTCGeometryTypeを渡します。
このenum値によって、この「3D空間上のオブジェクト」がどういう形式のバッファとして確保されるかをEmbree3のKernel側に伝えます。
例えばRTC_GEOMETRY_TYPE_TRIANGLEを渡せば、ポリゴンメッシュ形式(Vertex Buffer & Index Buffer)でメモリ上に確保されることを指定できます。

バッファの実体確保とジオメトリへのバインド

次に、「3D空間上のオブジェクト」の実体を確保してジオメトリにバインドします。
実体の動的確保をEmbree3の関数にやってもらう方法と、ユーザ自身で確保してEmbree3のKernel側に伝える方法があります。

順に説明します。

// embree3/rtcore_geometry.h
/*
  Embree3に動的確保をやってもらう方法
*/
void* rtcSetNewGeometryBuffer(
  RTCGeometry geometry,    // 上で生成したジオメトリ
  enum RTCBufferType type, // バッファ形式
  unsigned int slot,       // 
  enum RTCFormat format,   // 要素1つが持つ変数の型と数の指定
  size_t byteStride,       // 要素1つ分のサイズ
  size_t itemCount         // 要素数
);

struct Vertex {
  float x, y, z;
};
struct PolygonIndex{
  unsigned int v0, v1, v2;
}
// ~中略~
constexpr size_t numVertex = 6;
Vertex* geometryVertices = reinterpret_cast<Vertex*>(rtcSetNewGeometryBuffer(geometryHandle,
                                                                               RTC_BUFFER_TYPE_VERTEX,
                                                                               0,
                                                                               RTC_FORMAT_FLOAT3,
                                                                               sizeof(Vertex),
                                                                               numVertex));

Index* geometryIndices = reinterpret_cast<PolygonIndex*>(rtcSetNewGeometryBuffer(geometryHandle,
                                                                          RTC_BUFFER_TYPE_INDEX,
                                                                          0,
                                                                          RTC_FORMAT_UINT3,
                                                                          sizeof(PolygonIndex),
                                                                          numPolygons));

rtcSetNewGeometryBufferを利用すると、Embree3が引数からサイズを計算してバッファ領域の動的確保をしてくれます。
上記の例では、返り値のvoid *をキャストしてgeometryVerticesで受け取り、今後配列として利用します。

// embree3/rtcore_geometry.h
/*
  ユーザ自身で確保してEmbree3のKernel側に伝える方法
*/
void rtcSetSharedGeometryBuffer(
  RTCGeometry geometry,
  enum RTCBufferType type,
  unsigned int slot,
  enum RTCFormat format,
  const void* ptr,
  size_t byteOffset,
  size_t byteStride,
  size_t itemCount
);

struct PolygonVertex {
  float x, y, z;
};
struct PolygonIndex{
  unsigned int v0, v1, v2;
}
// ~中略~
constexpr size_t numVertex = 6;
constexpr int numPolygons = 8;
std::array<PolygonVertex, numVertex> geometryVertices;
std::array<PolygonIndex, numPolygons> geometryPolygons;

rtcSetSharedGeometryBuffer(geometryHandle,
                           RTC_BUFFER_TYPE_VERTEX,
                           0,
                           RTC_FORMAT_FLOAT3,
                           &geometryVertices,
                           0,
                           sizeof(Vertex),
                           geometryVertices.size());
rtcSetSharedGeometryBuffer(geometryHandle,
                           RTC_BUFFER_TYPE_INDEX,
                           0,
                           RTC_FORMAT_UINT3,
                           &geometryPolygons,
                           0,
                           sizeof(PolygonIndex),
                           geometryPolygons.size());

rtcSetSharedGeometryBufferを利用すると、既にユーザ側で確保したバッファ領域をジオメトリの実体として利用できます。
上記の例では、std::array<PolygonVertex, numVertex>として定義したgeometryVerticesのアドレスを引数に渡します。

バッファへのデータ入力

次に、上のいずれかの方法で定義したジオメトリに頂点データを入力します。
今回は手入力としますが、この処理を3Dモデルデータ読込処理に差し替えてしまえば、インポートが実現できます。
rtcSetSharedGeometryBufferを利用してバッファ領域確保を行なう場合は、geometryVerticesの定義時に以下のデータ入力処理をしておいても問題ありません。

// Embree3の関数は使用しません。

geometryVertices[0] = {-0.5f,  0.0f,  0.0f};
geometryVertices[1] = {-0.0f,  1.0f,  0.0f};
geometryVertices[2] = { 0.5f,  0.0f,  0.0f};
geometryVertices[3] = { 0.0f, -1.0f,  0.0f};
geometryVertices[4] = { 0.0f,  0.0f,  0.5f};
geometryVertices[5] = { 0.0f,  0.0f, -0.5f};
geometryPolygons[0] = { 1, 0, 4 };
geometryPolygons[1] = { 1, 4, 2 };
geometryPolygons[2] = { 3, 0, 4 };
geometryPolygons[3] = { 3, 4, 2 };
geometryPolygons[4] = { 1, 5, 0 };
geometryPolygons[5] = { 1, 2, 5 };
geometryPolygons[6] = { 3, 5, 0 };
geometryPolygons[7] = { 3, 2, 5 };

ジオメトリのコミット

ジオメトリのバッファへのデータ入力まで終えた後、このジオメトリの状態をEmbree3のKernel側へ伝えます。

// embree3/rtcore_geometry.h
void rtcCommitGeometry(RTCGeometry geometry);

rtcCommitGeometry(geometry);

今後geometryVertices(VertexBuffer)の各要素の値に少しでも変更があった場合、この関数を呼び出さなければその変更がEmbree3のKernel側に伝わらず、後のレイトレースの結果もコミット前と変わりません。
例えば、リアルタイムレンダリングで描画フレームごとににオブジェクトを移動させたいとき、毎フレームの更新処理でこの関数を呼び出す必要があります。

シーンへのジオメトリ登録

前ステップで初期化したジオメトリをシーンに登録します。
1つのシーンに複数のジオメトリを登録することが可能です。

// embree3/rtcore_scene.h
unsigned int rtcAttachGeometry(
  RTCScene scene,
  RTCGeometry geometry
);

rtcAttachGeometry(scene, geometry);

2019/12/14 typo修正 rtcAttachGeometry(geometry); から rtcAttachGeometry(scene, geometry);に修正しました。

登録の後、ジオメトリを解放することができます。

// embree3/rtcore_geometry.h
void rtcReleaseGeometry(RTCGeometry geometry);

rtcReleaseGeometry(geometry);

この関数呼び出しは、rtcAttachGeometryによってシーン、すなわちEmbree3のKernel側)がgeometryへの参照を獲得したため、ユーザが管理をする必要が無くなったから解放しようということです。
ただし、必ずしもこのタイミングではジオメトリ解放を行なう必要はありません。必要に応じてどうぞ。

シーン構築

シーンの状態をEmbree3のKernel側へ伝えます。

// declared in embree3/rtcore_scene.h
void rtcCommitScene(RTCScene scene);

rtcCommitScene(scene);

この関数が呼び出されると、Embree3は内部的にBVHの構築を行ないます。
ここまでで描画するシーンのセットアップがすべて完了となります。

レイトレース

1本のレイを飛ばす

// embree3/rtcore_common.h
struct RTCIntersectContext;
void rtcInitIntersectContext(RTCIntersectContext* context);

// embree3/rtcore_ray.h
struct RTCRay{
  float org_x;        // x coordinate of ray origin
  float org_y;        // y coordinate of ray origin
  float org_z;        // z coordinate of ray origin
  float tnear;        // start of ray segment

  float dir_x;        // x coordinate of ray direction
  float dir_y;        // y coordinate of ray direction
  float dir_z;        // z coordinate of ray direction
  float time;         // time of this ray for motion blur

  float tfar;         // end of ray segment (set to hit distance)
  unsigned int mask;  // ray mask
  unsigned int id;    // ray ID
  unsigned int flags; // ray flags
};
struct RTCHit
{
  float Ng_x;          // x coordinate of geometry normal
  float Ng_y;          // y coordinate of geometry normal
  float Ng_z;          // z coordinate of geometry normal

  float u;             // barycentric u coordinate of hit
  float v;             // barycentric v coordinate of hit

  unsigned int primID; // primitive ID
  unsigned int geomID; // geometry ID
  unsigned int instID[RTC_MAX_INSTANCE_LEVEL_COUNT]; // instance ID
};
struct RTCRayHit{
  struct RTCRay ray;
  struct RTCHit hit;
};
void rtcIntersect1(
    RTCScene scene, 
    RTCIntersectContext* context,
    RTCRayHit* rayhit
);

constexpr float FLOAT_INFINITY = std::numeric_limits<float>::max();

RTCIntersectContext context;
rtcInitIntersectContext(&context);

RTCRayHit rayhit;
rayhit.ray.org_x = 0.0f;
rayhit.ray.org_y = 0.0f;   // 始点(0, 0, -10)  
rayhit.ray.org_z = -10.0f; 
rayhit.ray.dir_x = 0.0f;
rayhit.ray.dir_y = 0.0f;   // 方向(0, 0, 1)
rayhit.ray.dir_z = 1.0f;
rayhit.ray.tnear = 0.0f;             // 始点から 
rayhit.ray.tfar = FLOAT_INFINITY; // 距離∞だけレイを飛ばす
rayhit.ray.flags = false; // このレイが一度何かと交差したか否か

rtcIntersect1(sceneHandle, &context, &rayhit);

if(rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID){
    // シーン上の何かしらと交差した
}
else{
    // 何にも交差しなかった
}

RTCRayHitのインスタンスを生成し、rtcIntersect1によって交差判定を行います。
すると、rayhit.hit.geomIDに、交差した最も近いジオメトリのIDが保存されます。
RTCRayHit::dir_xら方向ベクトル成分は正規化されている必要はありません。

スクリーンにレイを飛ばす

今回は簡単のため、カメラ位置などは考慮せず、平行投影と同様な交差判定をします。
100x100ピクセルのスクリーンを例として取り上げます。
screen screenxy

// 1本の場合と同様

constexpr float FLOAT_INFINITY = std::numeric_limits<float>::max();
constexpr uint32_t SCREEN_WIDTH_PX = 100;
constexpr uint32_t SCREEN_HEIGHT_PX = 100;
constexpr uint32_t HALF_HEIGHT = static_cast<int>(SCREEN_HEIGHT_PX) / 2;
constexpr uint32_t HALF_WIDTH = static_cast<int>(SCREEN_WIDTH_PX) / 2;

RTCIntersectContext context;
rtcInitIntersectContext(&context);

RTCRayHit rayhit;
rayhit.ray.org_z = -10.0f; // 始点は"手前"側
rayhit.ray.dir_x = 0.0f;   //
rayhit.ray.dir_y = 0.0f;   // (0, 0, 1)の方向へ飛ばす
rayhit.ray.dir_z = 1.0f;   //
rayhit.ray.tnear = 0.0f;

for (int32_t y = 0; y < SCREEN_HEIGHT_PX; ++y) {
  for (int32_t x = 0; x < SCREEN_WIDTH_PX; ++x) {
    rayhit.ray.org_x = (x - static_cast<int>(HALF_WIDTH)) / static_cast<float>(HALF_WIDTH);  // 毎ループ始点位置を更新
    rayhit.ray.org_y = (y - static_cast<int>(HALF_HEIGHT)) / static_cast<float>(HALF_HEIGHT); // 毎ループ始点位置を更新
    rayhit.ray.tfar = FLOAT_INFINITY;
    rayhit.ray.flags = false;
    rayhit.hit.geomID = RTC_INVALID_GEOMETRY_ID;
    rtcIntersect1(sceneHandle, &context, &rayhit);
    if (rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) {
      // シーン上の何かしらと交差した
    }
    else {
      // 何にも交差しなかった
    }
  }
}

ちょっとした注意点ですが、RTCRayHit::RTCRay::tfar, RTCRayHit::RTCRay::flags, RTCRayHit::RTCRay::geomIDはrtcIntersect1を呼び出すたびに書き込まれて変更される可能性があります。
ループ中で1つのRTCRayHitインスタンスを使いまわす場合は毎回のrtcIntersect1実行前にしっかり元の値に戻してあげましょう。

終了処理

ユーザ側で管理しているEmbree3オブジェクトを全て解放します。

// declared in embree3/rtcore_scene.h
void rtcReleaseScene(RTCScene scene);

// declared in rtcore_device.h
void rtcReleaseDevice(RTCDevice device);

rtcReleaseGeometry(geometryHandle); // 以前にやっていなければ
rtcReleaseScene(sceneHandle);
rtcReleaseDevice(deviceHandle);

まとめ

同様のシーンに対して平行投影的に交差判定し、100x100ピクセルのppm画像に書き込みする例です。
部分部分はこれまでのコードとほとんど同じです。

ソースコード

#include <embree3/rtcore.h>
#include <array>
#include <vector>
#include <limits>
#include <fstream>

struct Vertex {
  float x, y, z;
};

struct PolygonIndex {
  unsigned int v0, v1, v2;
};

struct ColorRGB {
  uint8_t r, g, b;
};

constexpr uint32_t SCREEN_WIDTH_PX = 100;
constexpr uint32_t SCREEN_HEIGHT_PX = 100;

void SaveAsPpm(const std::vector<std::array<ColorRGB, SCREEN_WIDTH_PX>>& image, const char* const fileName) {
  std::ofstream ofs(fileName);

  if (!ofs) {
    std::exit(1);
  }

  ofs << "P3\n"
      << SCREEN_WIDTH_PX << " " << SCREEN_HEIGHT_PX << "\n255\n";
  for (const auto& row : image) {
    for (const auto& pixel : row) {
      ofs << static_cast<int>(pixel.r) << " "
          << static_cast<int>(pixel.g) << " "
          << static_cast<int>(pixel.b) << "\n";
    }
  }

  ofs.close();
}

int main() {
  RTCDevice deviceHandle = rtcNewDevice(nullptr);
  RTCScene sceneHandle = rtcNewScene(deviceHandle);
  RTCGeometry geometryHandle = rtcNewGeometry(deviceHandle, RTC_GEOMETRY_TYPE_TRIANGLE);

  constexpr size_t numVertices = 6;
  constexpr size_t numPolygons = 8;

  std::array<Vertex, numVertices> geometryVertices;
  std::array<PolygonIndex, numPolygons> geometryPolygons;

  rtcSetSharedGeometryBuffer(geometryHandle,
                             RTC_BUFFER_TYPE_VERTEX,
                             0,
                             RTC_FORMAT_FLOAT3,
                             &geometryVertices,
                             0,
                             sizeof(Vertex),
                             geometryVertices.size());

  geometryVertices[0] = { -0.5f,  0.0f,  0.0f };
  geometryVertices[1] = {  0.0f,  1.0f,  0.0f };
  geometryVertices[2] = {  0.5f,  0.0f,  0.0f };
  geometryVertices[3] = {  0.0f, -1.0f,  0.0f };
  geometryVertices[4] = {  0.0f,  0.0f, -0.5f };
  geometryVertices[5] = {  0.0f,  0.0f,  0.5f };

  rtcSetSharedGeometryBuffer(geometryHandle,
                             RTC_BUFFER_TYPE_INDEX,
                             0,
                             RTC_FORMAT_UINT3,
                             &geometryPolygons,
                             0,
                             sizeof(PolygonIndex),
                             geometryPolygons.size());

  geometryPolygons[0] = { 1, 0, 4 };
  geometryPolygons[1] = { 1, 4, 2 };
  geometryPolygons[2] = { 3, 0, 4 };
  geometryPolygons[3] = { 3, 4, 2 };
  geometryPolygons[4] = { 1, 5, 0 };
  geometryPolygons[5] = { 1, 2, 5 };
  geometryPolygons[6] = { 3, 5, 0 };
  geometryPolygons[7] = { 3, 2, 5 };

  rtcCommitGeometry(geometryHandle);
  rtcAttachGeometry(sceneHandle, geometryHandle);
  rtcReleaseGeometry(geometryHandle);
  rtcCommitScene(sceneHandle);

  constexpr float FLOAT_INFINITY = std::numeric_limits<float>::max();
  std::vector<std::array<ColorRGB, SCREEN_WIDTH_PX>> screen{SCREEN_HEIGHT_PX};

  RTCIntersectContext context;
  rtcInitIntersectContext(&context);

  RTCRayHit rayhit;
  rayhit.ray.org_z = -10.0f;
  rayhit.ray.dir_x = 0.0f;
  rayhit.ray.dir_y = 0.0f;
  rayhit.ray.dir_z = 1.0f;
  rayhit.ray.tnear = 0.0f;

  constexpr uint32_t HALF_HEIGHT = static_cast<int>(SCREEN_HEIGHT_PX) / 2;
  constexpr uint32_t HALF_WIDTH = static_cast<int>(SCREEN_WIDTH_PX) / 2;

  for (int32_t y = 0; y < SCREEN_HEIGHT_PX; ++y) {
    for (int32_t x = 0; x < SCREEN_WIDTH_PX; ++x) {
      rayhit.ray.org_x = (x - static_cast<int>(HALF_WIDTH)) / static_cast<float>(HALF_WIDTH);
      rayhit.ray.org_y = (y - static_cast<int>(HALF_HEIGHT)) / static_cast<float>(HALF_HEIGHT);
      rayhit.ray.tfar = FLOAT_INFINITY;
      rayhit.ray.flags = false;
      rayhit.hit.geomID = RTC_INVALID_GEOMETRY_ID;
      rtcIntersect1(sceneHandle, &context, &rayhit);
      if (rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) {
        screen[y][x].r = 0xff;
        screen[y][x].g = 0x00;
        screen[y][x].b = 0xff;
      }
      else {
        screen[y][x].r = 0x00;
        screen[y][x].g = 0x00;
        screen[y][x].b = 0x00;
      }
    }
  }

  SaveAsPpm(screen, "result.ppm");

  rtcReleaseScene(sceneHandle);
  rtcReleaseDevice(deviceHandle);
}

結果

3Dモデルのシルエット表示

おわりに

お疲れさまでした。
次回解説でPhong Shadingの実装解説を予定していましたが体力が尽きました
もしこの記事内に間違いや不適切な記述を見つけた場合は、コメントでお教えください。

参考記事

EmbreeによるCPUレンダリング(その１) 交差判定
Embree公式ページ(Overview)
Embree公式ページ(API reference)
Embree公式リポジトリ
Wikipedia BVH

C++ SSE/AVX 入門の記録

※ この記事は、筆者が以前Qiitaに投稿したものの再掲版です。※

はじめに

Intelが考案したSSE/AVXという命令セットを利用することで、CPUで実行するコードのSIMD化によって高速化が図れます。
SIMD化はどういうものかから解説を始め、どうやってどのくらい高速化できたかなどを備忘録を兼ねて記していきます。
勉強中ゆえどうしても誤った記述をしてしまうかもしれません。見つけた方はコメントにて教えていただけると大変助かります！
また、この記事は Aizu Advent Calendar 2019 の5日目の記事です。
興味がある方は覗いてみてください。

想定読者

  • 基本的なC/C++の文法が分かる方
  • SIMDとかSSEって何？単語だけは耳にしたことあるけど...くらいの方。
  • ごく基本的なCPU周辺のアーキテクチャを知っている方。CPUは演算ユニット、レジスタ、キャッシュがあってRAMと併せて値を読み書きしながら実行して... 程度がわかっていればOK

概観の説明

SISDとSIMD

CPUは実行ファイル(や、アセンブリプログラム)に記述された命令を読み込んで実行します。
SSE/AVXを利用しない通常の命令セットでは、1つのデータに対して1つの演算を行います。つまりフリンの分類)で言うSISDというものです。 すると、CPUリソースが無駄になってしまうことがあります。
(2019/12/07追記) また任意の演算において、計算する必要のない領域が生まれることがあります。これではCPUリソースの無駄になってしまいます。

(注釈) 曖昧な書き方をしてしまっていたので訂正します。 「SISD演算が原因でCPUリソースが無駄になる」というわけではありません。

例えばこんな感じです。
演算で使うストレージ
この無駄な部分を利用して、本来次のステップで演算される予定だった別のデータに対する演算をこのステップで同時に実行できるようにする、というのがSIMD演算のアイデアです。

SSE/AVX

CPUでのSIMD演算を可能にしたx86命令セットの拡張が、Intelが考案したSSE/AVXです。
それぞれ Streaming SIMD Extensions, Advanced Vector Extensionsといいます。AVXは後に考案されたSSEの後継バージョンです。
これらはそれぞれ128bitsのxmmレジスタ、256bitsのymmレジスタを利用して演算します。
このレジスタはSIMD演算のためにCPUへ組み込まれた追加のレジスタで、それ意外の用途には使われません。
ymmレジスタはymm0~ymm15の16個ぶんあります。xmmレジスタも同数で、ymmの半分の領域を共有して利用します。
xmm/ymmレジスタ図解

C言語のIntrinsics

プログラマーがCPU命令セットを実行するよう記述するには普通アセンブリ言語でプログラムを書く必要がありますが、それをC言語から呼び出せるようにした組み込み関数が用意されています。
これらはIntrinsicsと呼ばれる関数群です。
SSE/AVXも同様にIntrinsicとしてC言語から実行できます。

tip : 並列処理とSIMD演算の違い

並列処理は、1つの処理をCPUの複数のスレッドで分担させて行います。
対して、SIMD演算は複数のデータに対する1つの命令をCPUの1つのコア内のxmm/ymmレジスタで分担させて行います。
マルチスレッド処理/SIMD演算図解

C++プログラムのSIMD化

C++によるIntrinsics実行

Intrinsicsを利用する際は、

  1. メモリからxmm/ymmレジスタへ値をロード
  2. レジスタ間演算
  3. レジスタからメモリへ値をストア

という流れで実行します。これらはそれぞれ

__m256 _mm256_load_ps (float const * mem_addr)
__m256 _mm256_add_ps (__m256 a, __m256 b)
void _mm256_store_ps (float * mem_addr, __m256 a)

などの関数を呼び出すことで実行できます。
xmm/ymmレジスタに格納する値はuint8 * 32, float * 8, double * 4などプログラマーが好きに切って使えます。
そのため、各型に対応するIntrinsicsが大量に用意されています。
Intel® Intrinsics Guideがそのリファレンスサイトとなっているので、関数を探しながら書くのがよさそうです。

例 : ベクトルと行列の積

お膳立てが済んだところで、C++からのIntrinsics呼び出しによる高速化を試してみます。
要旨と関係ない部分のコードは冗長になるので割愛します。ご容赦を...

(筆者が書いたプログラム全体を見たい方はこちらのリポジトリへ)

Vector3f Multiply(const Vector3f& vec, const Matrix4x4& matrix) noexcept {
  return {
    vec.x * matrix[0][0] + vec.y * matrix[1][0] + vec.z * matrix[2][0] + matrix[3][0], // 乗算3回 + 加算3回
    vec.x * matrix[0][1] + vec.y * matrix[1][1] + vec.z * matrix[2][1] + matrix[3][1], // 乗算3回 + 加算3回
    vec.x * matrix[0][2] + vec.y * matrix[1][2] + vec.z * matrix[2][2] + matrix[3][2], // 乗算3回 + 加算3回
  };
}

Vector3f MultiplySimd(const Vector3f& vec, const Matrix4x4& matrix) noexcept {
  __m256 row01Reg = _mm256_load_ps(reinterpret_cast<const float*>(&(matrix[0])));
  __m256 row23Reg = _mm256_load_ps(reinterpret_cast<const float*>(&(matrix[2])));
  __m256 multFactorFor01 = _mm256_set_ps(vec.x, vec.x, vec.x, vec.x, vec.y, vec.y, vec.y, vec.y);
  __m256 multFactorFor23 = _mm256_set_ps(vec.z, vec.z, vec.z, vec.z, 1.0f, 1.0f, 1.0f, 1.0f);
  row01Reg = _mm256_mul_ps(row01Reg, multFactorFor01);  // 乗算1回
  row23Reg = _mm256_mul_ps(row23Reg, multFactorFor23);  // 乗算1回
  __m256 resultReg = _mm256_add_ps(row01Reg, row23Reg); // 加算1回
  alignas(32) float resultMem[8];
  _mm256_store_ps(reinterpret_cast<float*>(resultMem), resultReg);
  return { resultMem[0] + resultMem[4], resultMem[1] + resultMem[5], resultMem[2] + resultMem[6] }; // 加算3回
}

ざっくりカウントすると

  • 普通に実装 : 乗算9回、加算9回
  • SIMD化実装 : 乗算2回、加算4回 + ymmレジスタへのロード&ストアコスト

演算回数が少なくなっているので、高速化が期待できます。

ベンチマーク

前項の内容を含む時間計測の実装です。

#include "randomgenerator.hpp"
#include "stopwatch.hpp"
#include "vector3.hpp"
#include "vector4.hpp"
#include "matrix4x4.hpp"

#include <functional>
#include <iostream>
#include <string>
#include <vector>

template <class ReturnTy, class LhsTy, class RhsTy>
void Benchmark(const std::vector<LhsTy>& lhs, const std::vector<RhsTy>& rhs, const size_t numCalc, std::function<ReturnTy(const LhsTy&, const RhsTy&)> Func, const char* const benchName) {
  // 普通のストップウォッチ(中身はただのstd::chrono関連関数呼び出し)
  StopWatch stopwatch;
  stopwatch.Start();
  for (auto i = 0u; i < numCalc; ++i) {
    Func(lhs[i], rhs[i]);
  }
  stopwatch.Stop();
  std::cout << benchName << " : " << stopwatch.GetMicroseconds() << " microsec" << std::endl;
  stopwatch.Reset();
}

int main(int argc, char* argv[]) {
  // ベンチマーク対象の関数を実行する回数
  const auto numCalc = (argc > 1) ? std::stoul(argv[1]) : 100000;
  std::cout << "<<" << numCalc << " times calculations as Benchmark>>" << std::endl;

  // [-100.0, 100.0]の範囲の一様分布に準ずる疑似乱数生成器
  RandomGenerator rand;

  // 時間計測 : 3次元ベクトルの内積・外積
  {
    std::vector<Vector3f> vecA(numCalc);
    std::vector<Vector3f> vecB(numCalc);
    for (auto i = 0u; i < numCalc; ++i) {
      vecA[i] = { rand(),
                  rand(),
                  rand() };
    }
    Benchmark<float, Vector3f, Vector3f>(vecA, vecB, numCalc, Dot3D, "normal Vector3 Dot");
    Benchmark<float, Vector3f, Vector3f>(vecA, vecB, numCalc, Dot3DSimd, "SIMD Vector3 Dot");
    Benchmark<Vector3f, Vector3f, Vector3f>(vecA, vecB, numCalc, Cross3D, "normal Vector3 Cross");
    Benchmark<Vector3f, Vector3f, Vector3f>(vecA, vecB, numCalc, Cross3DSimd, "SIMD Vector3 Cross");
  }

  // 時間計測 : 4次元ベクトルの内積
  {
    std::vector<Vector4f> vecA(numCalc);
    std::vector<Vector4f> vecB(numCalc);
    for (auto i = 0u; i < numCalc; ++i) {
      vecA[i] = { rand(),
                  rand(),
                  rand(),
                  rand() };
    }
    Benchmark<float, Vector4f, Vector4f>(vecA, vecB, numCalc, Dot4D, "normal Vector4 Dot");
    Benchmark<float, Vector4f, Vector4f>(vecA, vecB, numCalc, Dot4DSimd, "SIMD Vector4 Dot");
  }

  // 時間計測 : 3次元ベクトルと4次元行列の積
  {
    std::vector<Vector3f> vec(numCalc);
    std::vector<Matrix4x4> matrices(numCalc);
    for (auto i = 0u; i < numCalc; ++i) {
      vec[i] = { rand(),
                 rand(),
                 rand() };
      matrices[i] = { { rand(), rand(), rand(), rand() },
                      { rand(), rand(), rand(), rand() },
                      { rand(), rand(), rand(), rand() },
                      { rand(), rand(), rand(), rand() } };
    }
    Benchmark<Vector3f, Vector3f, Matrix4x4>(vec, matrices, numCalc, Multiply, "normal Vector*Matrix");
    Benchmark<Vector3f, Vector3f, Matrix4x4>(vec, matrices, numCalc, MultiplySimd, "SIMD Vector*Matrix");
  }
}

実行結果:

g++ -O3 オプション

<<100000 times calculations as Benchmark>>
normal Vector3 Dot : 232 microsec
SIMD Vector3 Dot : 245 microsec
normal Vector3 Cross : 797 microsec
SIMD Vector3 Cross : 850 microsec
normal Vector4 Dot : 223 microsec
SIMD Vector4 Dot : 250 microsec
normal Vector*Matrix : 2179 microsec
SIMD Vector*Matrix : 1162 microsec

g++ -O0 オプション

<<100000 times calculations as Benchmark>>
normal Vector3 Dot : 2412 microsec
SIMD Vector3 Dot : 3280 microsec
normal Vector3 Cross : 4755 microsec
SIMD Vector3 Cross : 6964 microsec
normal Vector4 Dot : 2486 microsec
SIMD Vector4 Dot : 3172 microsec
normal Vector*Matrix : 9324 microsec
SIMD Vector*Matrix : 7598 microsec

3次元ベクトルの内積/外積程度の演算回数ではあまり高速化できないどころか、(おそらく)ロード&ストアのオーバーヘッドのせいで遅くなってしまうようです。
対して、演算回数の多いベクトルと行列の積では-O3で2倍、-O0で1.2倍程度高速化できているのがわかります。

実行環境など

CPU : AMD Ryzen 2700X
OS  : 18.04.3 LTS (Bionic Beaver)
compiler : gcc 7.4.0

注意点・ハマったところなど

データのアラインメント

lord元/store先になるメモリ領域がSSEなら16bytes-aligned, AVXなら32bytes-alignedになっていないとSegmentation Faultでプログラムが落ちる。

  __m128 left = _mm_set_ps(0.0f, 1.0f, 2.0f, 3.0f);
  __m128 right = _mm_set_ps(0.0f, 1.0f, 2.0f, 3.0f);
  __m128 resultReg = _mm_mul_ps(left, right);
  alignas(16) float resultMem[4];  // alignas(16)が無いと環境によってはセグる!!
  _mm_store_ps(reinterpret_cast<float*>(resultMem), resultReg);

分からないところ・疑問など

一通り触ってみて、記事執筆中に湧いた新たな疑問などを参考のため残します。
これらを解説する資料や技術書など,もしあればぜひコメントで教えてください！

  • xmm/ymmレジスタはCPUの1コアにつき1セット(=16個のレジスタ)存在している？
  • もし ↑ の通りならば、Intrinsicsを含むコードを実行する並列スレッド数(std::threadなどのインスタンス数)をCPUのコア数より多くすればxmm/ymmレジスタの値が破綻し得る？
  • コンパイラの最適化は、可能な部分を自動でSIMD化してくれるのか？
  • SSE/AVXのほかにFMAという命令セット拡張があり、SSE->AVXのようにサイズが増えたわけではないらしい。用途は何だろう？
  • SSE/AVX Intrinsics利用時はプログラマーが手動でレジスタ-メモリ間のロード/ストアをするが、普段C++で普通に書く int n = a + b; のような演算もレジスタを利用して演算するはず。それらはどのような流れでレジスタに到達し、CPUのコアの演算ユニットで演算され、再びメモリに戻ってくるのか？
  • 普通の演算、SIMD化演算問わずキャッシュはどのようにプログラム実行に関係してくるのか？

おわりに(感想)

少し前からx86アセンブリ言語を勉強すべきか？とは思っていましたが、SSE/AVX Intrinsicsに触れてみてその機運がとても高まりました。
ここまで読んでいただいた方、ありがとうございました。何かの参考になれば幸いです。

参考文献

Intel AVX を使用して SIMD 演算を試してみる
Intel AVX を使用して SIMD 演算を試してみる - その2 -
Introduction to Intel® Advanced Vector Extensions
Intel® Intrinsics Guide
組み込み関数（intrinsic）によるSIMD入門
SIMD - primitive: blog

GitHub Actions 導入までの流れ

GitHubで利用できるCI/CD機能 GitHub Actions を利用する際の流れまとめ。
備忘録的に記していきます。

1. YAMLファイル記述と設置

Actionsを利用するために、YAMLファイルに自動実行するプロセスを記述します。

YAMLファイルの置き場所

リポジトリのルートディレクトリ直下に.github/workflowsというディレクトリを作り、その下に任意のファイル名でYAMLファイルを設置します。
すると、Actionsが自動で読み込んで特定のタイミングで実行してくれるようになります。

# 例
REPOSITORY_ROOT/.github/workflows/build_verification.yml

記述項目

Actions固有の設定項目や実行プロセス名などを書き込んでいきます。
詳しくは -> GitHub Actionsのワークフロー構文 - GitHubヘルプ
以下重要なものを列挙します。

  • name : このYAMLファイルで実行するワークフローの名前指定。

  • on : ワークフローを実行するトリガーイベントの指定。

  • jobs : 一連のワークフロー実行のブロック。OSごとのテストや他のプロセスのひとまとまりなどをまとめる。この下(インデント的に)に↓の構文を記述していく。

    • <任意のジョブ名> : ジョブ名指定。
      • runs-on : ワークフローの各ジョブを実行する仮想環境の指定。仮想マシンについては後述。
      • needs : このジョブを実行する前に終わらせる必要のある他ジョブを指定。
      • steps : 1つ以上のアクションまたはシェルコマンドを実行。複数の場合は配列として記述する。
        • name : 1つのステップの名前。GitHub上で表示される。
        • uses : ステップの一部として実行するアクションを指定。公開されているパブリックリポジトリのものなどを使用できる。
        • run : シェルコマンドを実行する。
        • shell : 実行するシェルを指定。bashやpythonなど。bashはWin/mac/Linuxすべてで利用可。
      • container : ジョブを実行するコンテナの指定。
        • image : Dockerイメージの指定。Docker Hubイメージ名かpublic dockerレジストリ名で。
  # 例
  name : Build
  on   : [push, pull_request]
  jobs : 
    example_job:
      runs-on : ubuntu-latest
      
    release_build:
      runs-on : ubuntu-latest
      needs   : example_job
      steps   : 
        - name : Print path
          run  : echo $PATH
          shell: bash
          
        - uses : actions/setup-node@74bc508

2. ステータスバッジ設置

ワークフローのステータスを表示するバッジを表示します。
README.mdによくある Build|passingみたいなやつ。
以下のうちどちらかのURLで画像を取得できます。

https://github.com/<OWNER>/<REPOSITORY>/workflows/<WORKFLOW_NAME>/badge.svg
https://github.com/<OWNER>/<REPOSITORY>/workflows/<WORKFLOW_FILE_PATH>/badge.svg

あとはREADME.mdなどに
![build status](<上のURL>)
などとすればOK。

3. GitHub上で確認

私のリポジトリにこのようなYAMLファイルを追加し、確認してみます。

name: build

on: 
  push:
    branches:
      - autobuild_setup

jobs:
  build:
    runs-on: ubuntu-latest
    steps  :
      - name : test print
        run  : echo "this is test print"
        shell: bash

YAMLファイルをpushしたリポジトリのActionsタブを見ると、実行されたワークフローの履歴が表示されています。 github tab 最新のbuildを見てみると、こうなりました。
しっかりbuildのtest printが実行されてthis is test printが出力されています。

build detail

ステータス表示はこんな感じで。

github badge

appendix 1. ワークフロー実行の仮想環境

ワークフローを実行する仮想マシン環境を指定できます。

GitHubホストランナー

GitHubが標準で用意している環境を「GitHubホストランナー」といい、Actionsを実行するアプリケーションがインストールされています。
Windows macOS Linux のActions環境を用意してくれているようです。
LinuxとmacOSに関しては、パスワード無しでsudoつきコマンドを実行できます。
以下、yamlのruns-onに指定できるラベル一覧。

  • windows-latest : Windows Server 2019
  • ubuntu-latest : Ubuntu 18.04
  • ubuntu-18.04 : Ubuntu 18.04
  • macos-latest : macOS Catalina 10.15

これらの環境に標準でインストールされているコマンド/ツールはここ -> GitHubホストランナーにインストールされるソフトウェア - GitHubヘルプ

自己ホストランナー

実行する仮想マシン指定は、GitHubホストランナーでないものも利用できます。
この記事では詳しくは触れません。
以下参考リンク↓
自己ホストランナーについて - GitHubヘルプ
ワークフローでのセルフホストランナーの利用 - GitHubヘルプ

appendix 2. Dockerコンテナ

CI/CDを実現したいリポジトリとは別にリポジトリを作成し、そちらのルートディレクトリにDockerfile, action.yml, extrypoint.shを置くことで1つのアクションとして公開できます。
こうするとCI/CDを実現したいリポジトリのワークフローYAMLファイルの jobs.<jobname>.steps.usesで指定できるようになります。
詳細 -> Dockerコンテナのアクションを作成する - GitHubヘルプ

参考文献

ワークフローを設定する - GitHubヘルプ
GitHubホストランナーの仮想環境 - GitHubヘルプ
GitHub Actionsのワークフロー構文 - GitHubヘルプ
Dockerコンテナのアクションを作成する - GitHubヘルプ
GitHubホストランナーにインストールされるソフトウェア - GitHubヘルプ
自己ホストランナーについて - GitHubヘルプ
ワークフローでのセルフホストランナーの利用 - GitHubヘルプ

2020/1/1 追記
詳しい解説記事を見つけましたので追加しておきます。