粒子のレンダリング (2) ポイントの移動

グチが多い

このブログは「グチが多い」そうなんですけど、自分としては常にネタを仕込むことを心掛けております。しかし、久しぶりに書いた前回の記事はのっけから愚痴になってしまっていて、「うむむこれはまずい、やはり表向きにはリア充のふりをしておかなければならぬ」と思い直しております。あっ、でもリア充とかネアカとかネクラとか最近聞かんぞ、もしかしたら既に死語なのか？あっそうだ、関係ないけど（いやあるけど）この一連の記事は最初ぐっちに読まそうと思って書き始めたんだからね (グチだけに)。ちゃんと読めよ。

コンピュートシェーダ

ぐっちは粒子の運動に関して CPU でやるみたいな前提だったけど、やっぱりシェーダを使わんと色々限界があると思うんですよ。そこでコンピュートシェーダですよ（たかしはコンピュートシェーダでやろうとして四苦八苦してるみたいで、たった今も「どうやってデバッグしたらいいですか」とか聞きに来たけど）。なので、とにかくコンピュートシェーダをこのサンプルプログラムに組み込んでみます。まず、Blob クラスにコンピュートシェーダのプログラムオブジェクト名を保持するメンバと、粒子の位置を更新するメソッドを追加します。Blob.h を次のように書き換えます。

//
// 粒子群オブジェクト
//
class Blob
{
  // 頂点配列オブジェクト名
  GLuint vao;
 
  // 頂点バッファオブジェクト名
  GLuint vbo;
 
  // 頂点の数
  const GLsizei count;
 
  // 描画用のシェーダ
  const GLuint drawShader;
 
  // uniform 変数の場所
  const GLint mpLoc, mvLoc;
  
  // 更新用のシェーダ
  const GLuint updateShader;
  
public:
  
  ...
  
  // 描画
  void draw(const GgMatrix &mp, const GgMatrix &mv) const;
  
  // 更新
  void update() const;
};

次に、Blob.cpp を以下のように変更します。コンストラクタで追加したメンバに初期値を設定します。シェーダのソースプログラム update.comp については後述します。ggLoadComputeShader() 関数は授業の宿題の補助プログラムで用意しています。詳細はドキュメント(PDF) を参照してください。

// コンストラクタ
Blob::Blob(const Particles &particles)
  : count(static_cast<GLsizei>(particles.size()))
  , drawShader(ggLoadShader("point.vert", "point.frag"))
  , mpLoc(glGetUniformLocation(drawShader, "mp"))
  , mvLoc(glGetUniformLocation(drawShader, "mv"))
  , updateShader(ggLoadComputeShader("update.comp"))
{
  ...
}

シェーダストレージバッファオブジェクト

update() メソッドの実装は次のようにします。実際の処理はコンピュートシェーダで行います。描画に使うシェーダは glDrawArrays() や glDrawElements() などで描画を実行（ドローコール）する際に起動されますが、コンピュートシェーダは描画処理とは関係なく独立して起動することができます。コンピュートシェーダの起動には glDispatchCompute() を使います。

コンピュートシェーダの起動に先立って、コンピュートシェーダの入出力データを準備します。ここでは、これにシェーダストレージバッファオブジェクト (Shader Storage Buffer Object, SSBO) を用います。これは頂点バッファオブジェクト（ここでは vbo）を glBindBufferBase() により GL_SHADER_STORAGE_BUFFER に結合したものです。glBindBufferBase() の第 2 引数は 結合ポイント (Binding Point, BP) と呼ばれ、コンピュートシェーダ側でこの番号を指定します。ここでは 0 番を指定しています。

// 描画
void Blob::draw() const
{
  // 描画する頂点配列オブジェクトを指定する
  glBindVertexArray(vao);
  
  // 点で描画する
  glDrawArrays(GL_POINTS, 0, count);
}
 
// 更新
void Blob::update() const
{
  // シェーダストレージバッファオブジェクトを 0 番の結合ポイントに結合する
  glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, vbo);
  
  // 更新用のシェーダプログラムの使用開始
  glUseProgram(updateShader);
  
  // 計算を実行する
  glDispatchCompute(count, 1, 1);
}

これを粒子群オブジェクトの描画のあとに実行します。flow.cpp を以下のように変更します。

  ...
  
  //
  // 描画
  //
  while (window)
  {
    ...
 
    // 粒子群オブジェクトを描画する
    blob->draw();
 
    // 粒子群オブジェクトを更新する
    blob->update();
 
    // カラーバッファを入れ替える
    window.swapBuffers();
  }
}

コンピュートシェーダの起動

GPU の性能の最大のキモは、その並列処理の能力にあります。コンピュートシェーダによる処理も、同じ処理を行う複数のスレッドによって並列に実行されます。一つのスレッドが一つの CPU コアによる処理に相当します。ワークグループはこのスレッドを複数まとめたものであり、同じワークグループに属するスレッドは、シェアードメモリを介してデータを共有することができます。一つのワークグループに属するスレッドの数は、コンピュートシェーダのソースプログラムにおいて layout 修飾子により指定します。コンピュートシェーダのソースプログラムについては、後で説明します。

コンピュートシェーダのワークグループとスレッド

glDispatchCompute() の三つの引数 num_groups_x, num_groups_y, num_groups_z は、このワークグループを起動する数を指定します。これらはいずれも、少なくとも 65536 まで指定できます。実際に指定可能な数は、glGetIntegeri_v() の第 1 引数 pname に GL_MAX_COMPUTE_WORK_GROUP_COUNT を指定して調べることができます。第 2 引数が 0 なら num_groups_x、1 なら num_groups_y、2 なら num_groups_z に指定できる最大値が得られます。

NVIDIA GeForce GTX 1080Ti では num_groups_x = 2147483647 = $2^{31} - 1$、num_groups_y = num_groups_z = 65536 でした。また Intel HD Graphics 5000 (Core i5 4250U 内蔵) はいずれも 65536 でした。

例えば $640 \times 480$ 画素の画像を扱う場合、ワークグループのサイズを $32 \times 32 \times 1$ スレッドで構成すれば、glDispatchCompute() では $(640 / 32) \times (480 / 32) \times (1 / 1) = 20 \times 15 \times 1$ のワークグループを起動すればよいことになります。また、この後に示すコンピュートシェーダのソースプログラムでは、一つのワークグループあたりスレッドを一つだけ起動している (ワークグループのサイズが $1 \times 1 \times 1$) ため、粒子の数、すなわち count 個のワークグループを起動すればよいことになります。そのため update() メソッドでは、count $\times 1 \times 1$ のワークグループを起動しています。

Intel の HD Graphics 5000 のように num_groups_x が規格の最低値の 65536 しかないと、このようなプログラムでは割と簡単にこの制限に引っかかってしまいます。その場合は num_groups_x * num_groups_y * num_groups_z が count を超えるように設定し、コンピュートシェーダ内で GLSL の組み込み変数 gl_LocalInvocationIndex が count 未満の場合のみ計算を行うようにするなどの工夫が必要になります。

コンピュートシェーダのソースプログラム

この glDispatchCompute() によって起動するコンピュートシェーダのソースプログラムは、次のようなものです。ソースファイル名は update.comp とします。この最初の layout 修飾子の local_size_x、 local_size_y、 local_size_z には、同時に実行するシェーダプログラムの数を指定します。local_size_x * local_size_y * local_size_z 個のスレッドが（本当に並列に実行されるかどうかは別にして）一つのワークグループとして一斉に起動します。このサンプルプログラムでは、一つのワークグループあたり $1 \times 1 \times 1 = 1$ 個のスレッドが実行されることになります。またワークグループ自体も並列に動作しますが、実際にどれだけの数のワークグループが並列に動作するかどうかも、GPU の能力次第のようです。

また、粒子データのデータ型として、Particle.h で定義している Particle 構造体と同じ構造の GLSL の構造体を Particle という名前で宣言しておきます。そして、それを要素としたバッファ Particles を定義します。std430 はこのバッファのメモリレイアウトが C/C++ 言語に準じたものであることを示し、binding には参照するシェーダストレージバッファオブジェクトが結合されている結合ポイント (BP) を指定します。このサンプルプログラムでは 0 番に粒子データを格納したシェーダストレージバッファオブジェクトを結合していました。

一つのスレッドが参照する粒子データの番号は、このサンプルプログラムでは 1 ワークグループあたり 1 スレッドにしていたので、ワークグループの位置がそのままスレッドの位置になります。また update() メソッドでは x 方向に count 個のワークグループを起動してましたから、ワークグループの位置の x 座標がスレッドの番号になります。シェーダストレージバッファオブジェクトからその番号のデータを取り出して、データを更新します。ここでは粒子の位置の z 座標を 0.1 だけずらします。

#version 430 core
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;
 
// 粒子データ
struct Particle
{
  vec4 position;
  vec4 velocity;
};
 
// 粒子群データ
layout(std430, binding = 0) buffer Particles
{
  Particle particle[];
};
 
void main()
{
  // ワークグループ ID をそのまま頂点データのインデックスに使う
  const uint i = gl_WorkGroupID.x;
 
  // 位置を更新する
  particle[i].position.z += 0.1;
}

この local_size_x、local_size_y、local_size_z については、それぞれ少なくとも 1024, 1024, 64 の数が指定できます。この最大値も同様に GL_MAX_COMPUTE_WORK_GROUP_SIZE によって調べることができます。また、これらの三つの値の積、すなわちワークグループ内で起動可能なスレッドの数は、少なくとも 1024 あります。この最大数は glGetIntegeriv() の第 1 引数 pname に GL_MAX_COMPUTE_WORK_GROUP_INVOCATIONS を指定して調べることができます。

NVIDIA GeForce GTX 1080Ti では 1536 でした。

これらに加えて、コンピュートシェーダ内のシェアードメモリの合計サイズにも制限があります。これは少なくとも 32KB あります。実際に使用できるシェアードメモリのサイズは、同様に GL_MAX_COMPUTE_SHARED_MEMORY_SIZE によって調べることができます。

glDispatchCompute() の引数に指定した num_groups_x、num_groups_y、num_groups_z は、GLSL の組み込み変数 gl_NumWorkGroups に格納されています。このようにワークグループやスレッドが 3 次元の格子状に配置されているのは、実際に 3 次元的に処理を行うというより、一つのワークグループやスレッドが、自分の処理するデータの領域や個々の要素を把握できるようにすることが目的のようです。一つのワークグループのデータ全体における 3 次元の位置は、GLSL の組み込み変数 gl_WorkGroupID で調べることができます。また、一つのスレッドのワークグループ内での 3 次元の位置は、gl_LocalInvocationID で調べることができます。一つのワークグループのサイズ local_size_x、local_size_y、local_size_z は gl_WorkGroupSize で得られるので、一つのスレッドのデータ全体における位置は gl_WorkGroupID * gl_WorkGroupSize + gl_LocalInvocationID になります。この値は gl_GlobalInvocationID に格納されています。さらに、このワークグループ内を一次元に展開したときのスレッドの位置が gl_LocalInvocationIndex に入っています。

実行結果

Windows 10 ってウィンドウ開くときこういうアニメーションしてたのね…

コンピュートシェーダを追加したサンプルプログラム

重力をかけてみる

この粒子群に重力をかけてみます。粒子の位置を更新するシェーダプログラム update.comp を次のように修正します。重力の加速度ベクトルを gravity という uniform 変数に入れておきます。また、更新するタイムステップも dt という uniform 変数に入れておきます。なお、これらを uniform 変数にしているのは、あとでこれをシェーダプログラム外から変更するかもしれないからです。しないかもしれないです。

#version 430 core
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;
 
// 粒子データ
struct Particle
{
  vec4 position;
  vec4 velocity;
};
 
// 粒子群データ
layout(std430, binding = 0) buffer Particles
{
  Particle particle[];
};
 
// 重力加速度
uniform vec4 gravity = vec4(0.0, -9.8, 0.0, 0.0);
 
// タイムステップ
uniform float dt = 0.01666667;

粒子の現在の速度にタイムステップをかけたものを粒子の現在位置に加えて、粒子の位置を更新します。そのあと、重力加速度にタイムステップをかけたものを粒子の現在の速度に加えて、粒子の速度も更新します。これはオイラー法っていうやつです。

void main()
{
  // ワークグループ ID をそのまま頂点データのインデックスに使う
  const uint i = gl_WorkGroupID.x;
 
  // 位置を更新する
  particle[i].position += particle[i].velocity * dt;
 
  // 速度を更新する
  particle[i].velocity += gravity * dt;
}

地面で跳ね返らせてみる

重力をかけると下に落ちて行って見えなくなってしまうので、地面を付けようと思います。地面の高さと跳ね返ったときの速度の減衰率を、それぞれ uniform 変数 height と attenuation で設定し、もし地面に落ちたら跳ね返るようにしてみます。

#version 430 core
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;
 
...
 
// 重力加速度
uniform vec4 gravity = vec4(0.0, -9.8, 0.0, 0.0);
 
// 地面の高さ
uniform float height = -1.5;
 
// 減衰率
uniform float attenuation = 0.7;
 
// タイムステップ
uniform float dt = 0.01666667;
 
void main()
{
  // ワークグループ ID をそのまま頂点データのインデックスに使う
  const uint i = gl_WorkGroupID.x;
 
  // 位置を更新する
  particle[i].position += particle[i].velocity * dt;
 
  // 速度を更新する
  particle[i].velocity += gravity * dt;
 
  // もし地面に落ちたら
  if (particle[i].position.y < height)
  {
    // y 方向の速度を反転する
    particle[i].velocity.y = -attenuation * particle[i].velocity.y;

地面で跳ね返るので y 軸方向の速度の符号を反転し、それに減衰率 attenuation をかけます。ただし、このままだと地面に粒子がめり込んだまま、跳ね返っても二度と地表に出られないということになってしまいます。粒子の位置の更新は連続的に行っているわけではなく、画面の再表示タイミングに合わせて周期的に行っています。そのため、粒子が地面に落ちたと判定されたときの粒子の位置は、既に地面の下にあります。そのため、そこから粒子の向きを反転しても、減衰によってさらにその次のタイミングでも地面の下にあるということが起こります。

したがって、本当は粒子が地面で跳ね返ったあとの、次のタイミングでの位置を求める必要があります。しかし、めんどくさいので、粒子の高さ (y 座標値) を無理やり地面の高さにするという非常に安直な方法を採用します。こういうことをすると、本来は高さの異なる粒子が同じ高さにそろえられてしまって不自然に見えるため、良い子はここでちゃんと粒子の位置を計算するようにしましょう。宿題だ宿題！

    // 高さを地面の高さに戻す
    particle[i].position.y = height;
  }
}

ただですね、ここで動かしてみるとわかると思いますが、地面に落ちた粒子がいつまでもブルブル震えていると思います。粒子が地面付近で地表と地下を行ったり来たりしているのですね。そこで、ここでもう一つ安直な手段を採用します。速度を更新している部分と位置を更新している部分を入れ替え、速度を更新してから、その位置を使って位置を更新します。次のタイミングでの速度を使って現在の位置を更新するわけです。一応、これはシンプレクティック (Symplectic) 法っていう根拠のある方法らしいです。

void main()
{
  // ワークグループ ID をそのまま頂点データのインデックスに使う
  const uint i = gl_WorkGroupID.x;
 
  // 速度を更新する
  particle[i].velocity += gravity * dt;
 
  // 位置を更新する
  particle[i].position += particle[i].velocity * dt;
 
  // もし地面に落ちたら
  if (particle[i].position.y < height)
  {
    // 高さを地面の高さに戻す
    particle[i].position.y = height;
 
    // y 方向の速度を反転する
    particle[i].velocity.y = -attenuation * particle[i].velocity.y;
  }
}

これで落ち着くと思います。

繰り返し落とす

単位は繰り返し落としたりしないよう万全の注意を払っていただきたいものですが、このサンプルプログラムでは粒子は一度しか落とされないので、粒子が地面に落ちたあと止まってしまって面白くありません。そこで、粒子を繰り返し落とすようにしてみます。flow.cpp において定期的に粒子の初期値を設定する initialize() メソッドを呼び出すようにします。まず、繰り返しの間隔を決めます。

//
// アプリケーション本体
//
#include "GgApplication.h"
 
...
 
// 一つの粒子群の中心からの距離の標準偏差
const GLfloat pDeviation(0.3f);
 
// アニメーションの繰り返し間隔
const double interval(5.0);
 
...

そして図形の描画のループで経過時刻を調べ、それが繰り返しの間隔を超えていたら粒子の位置を元に戻し、時計をリセットします。関係ないけど、GLFW には今のところタイマーで関数を呼び出す機能はないみたいです。

  ...
  
  //
  // 描画
  //
  while (window)
  {
    // 定期的に粒子群オブジェクトをリセットする
    if (glfwGetTime() > interval)
    {
      blob->initialize(initial);
      glfwSetTime(0.0);
    }
 
    ...
 
    // カラーバッファを入れ替える
    window.swapBuffers();
  }
}

初期値として位置の代わりに速度を設定してみる

最後に、粒子群の生成の際には速度を設定するようにしてみます。粒子が多少派手に散らばるように、一つの粒子群の中心からの距離の標準偏差 pDeviation も大きくしてみます。

//
// アプリケーション本体
//
#include "GgApplication.h"
 
...
 
// 一つの粒子群の中心からの距離の標準偏差
const GLfloat pDeviation(1.0f);
 
// アニメーションの繰り返し間隔
const double interval(5.0);

粒子の位置は粒子群の中心位置 (cx, cy, cz) とし、もともと粒子群の中心からの相対位置として求めていたもの (r * sp * ct, r * sp * st, r * cp) を速度に用います。

//
// 粒子群の生成
//
//   particles 粒子群の格納先
//   count 粒子群の粒子数
//   cx, cy, cz 粒子群の中心位置
//   rn メルセンヌツイスタ法による乱数
//   mean 粒子の粒子群の中心からの距離の平均値
//   deviation 粒子の粒子群の中心からの距離の標準偏差
//
void generateParticles(Particles &particles, int count,
  GLfloat cx, GLfloat cy, GLfloat cz,
  std::mt19937 &rn, GLfloat mean, GLfloat deviation)
{
  ...
 
  // 原点中心に直径方向に正規分布する粒子群を発生する
  for (int i = 0; i < count; ++i)
  {
    // 緯度方向
    const GLfloat cp(uniform(rn) - 1.0f);
    const GLfloat sp(sqrt(1.0f - cp * cp));
 
    // 経度方向
    const GLfloat t(3.1415927f * uniform(rn));
    const GLfloat ct(cos(t)), st(sin(t));
 
    // 粒子の粒子群の中心からの距離 (半径)
    const GLfloat r(normal(rn));
 
    // 粒子を追加する
    particles.emplace_back(cx, cy, cz, r * sp * ct, r * sp * st, r * cp);
  }
}

こんな感じになります。これも一つの粒子群の中心からの距離の平均 pMean や一つの粒子群の中心からの距離の標準偏差 pDeviation を色々変えてみてください。

破裂して地面に落ちる粒子群

床井研究室

■ 2018年10月18日 - 粒子のレンダリング (2) ポイントの移動