水の表現の挑戦 multi threadで波の計算を高速化した話

描画される波形は前回から変わっていないが、ずっと気になっていたmulti threadで波のポリゴン(頂点のy座標、高さ)の計算を高速化してみた。

std::thread

使うライブラリはstd::threadです。
#include <thread>
でインクルードできます。

使い方

詳しくは
thread - cpprefjp C++日本語リファレンス

【C++】std::thread についての覚え書き - Qiita
をご覧ください。私はあるクラスの中でメンバ関数として呼び出して使用しているので

const uint32_t threadNum = 24;
std::array<std::unique_ptr<std::thread>, threadNum> threads;
for(uint32_t i = 0; i < threadNum; i++)
    threads[i] = std::make_unique<std::thread>(&Class::CalcSeaLevel, this, i, time, threadNum);
for(uint32_t i = 0; i < threadNum; i++)
    threads[i]->join();
for(uint32_t i = 0; i < threadNum; i++)
    threads[i].reset();
void Class::CalcSeaLevel(uint32_t mod, float time, uint32_t threadNum)
{
  uint32_t vertexSize = GetVertexSize();
  for(uint32_t i = 0 + mod; i < vertexSize; i = i + threadNum)
  {
      glm::vec3& pos = mVertices[i].pos;
      if(pos.y > -0.0001f)
          continue;
      pos.y = -1.5 + sin(time * 5.0f - glm::length(pos - glm::vec3(0, pos.y, 0)) * 8.0f) / 16.0f;
  }
}

としています。
上記の2つともクラスのメンバ関数ですが、1つ目の処理で2つ目の関数を呼び出しています。メンバ関数内からメンバ関数をポインタにして呼ぶ場合には

std::thread th(&Class::func, this, argument1, argument2, ...);

となることに注意。
c++ - Start thread with member function - Stack Overflow

処理してほしい関数を渡すだけのイメージなので、シンプルでよいですね。

この処理の後、15fps程度であったのが60fps程度にまで上がっていたので、処理はだいぶ高速化されるようです。

Simple Water Caustics



Vulkan Ray TracingでWater Causticsを書いてみました。Water Causticsとは水面の底にできる光の模様のことです。正確には光りの屈折でできあがる様々な模様のことのようです。
コースティクス - Wikipedia

水面の波の式

様々な波を試してみましたが、最も綺麗な模様を描けるのは円形波のようでした。
円形波自体は



y = A sin2\pi(\frac{t}{T}-\frac{r}{\lambda})
で表されます。ここでT=周期\lambda=波長r=\sqrt{x^2+z^2}です。

上記のプログラムでは

pos.y = -1.5 + sin(time * 5.0f - glm::length(pos - glm::vec3(0, pos.y, 0)) * 8.0f) / 16.0f;

としています。

Caustics描画

イデアは次のNVIDIAのblogに書かれているものと同じです。
Chapter 2. Rendering Water Caustics | NVIDIA Developer

f:id:Aqoole_Hateena:20210617231213p:plain
出典 https://developer.nvidia.com/gpugems/gpugems/part-i-natural-effects/chapter-2-rendering-water-caustics figure 2-8 Accessing the Environment Map

ただし今回の例だと光源が一つだけなのでlight mapは作成せずに、shaderの中で直接計算しています。

//caustics
uint causticFlags = gl_RayFlagsOpaqueEXT | gl_RayFlagsTerminateOnFirstHitEXT;
traceRayEXT(topLevelAS, causticFlags, 0xFF, 1, 0, 1, refractPos, 0.001, vec3(0.0, -1.0, 0.0), tMax, 2);
vec3 waterSurface = prd.pos;
vec3 waterNormal = prd.normal;
vec3 waterV = refract(vec3(0.0, -1.0, 0.0), -waterNormal, nWater/nAir);
float t = (pushC.lightPosition.y - waterSurface.y) / waterV.y;
vec3 sunMapPoint = waterSurface + t * waterV;
float caustics = 1.0 / (1.0 + length(sunMapPoint - pushC.lightPosition));
color += caustics;

図でいうところの



v1 = t * waterV
ということになります。


本来ならCausticsは光源から波に当たった光の内、屈折した結果ある点に集まるすべての光をトレースする必要がありますが、それだと難易度が高すぎるので、一番強い光ということで水面に垂直にあたる光のみをトレースしています。
これだけでも、実施例のようにそれなりの結果が得られますね。

水の表現の挑戦 波の表現について

上からの図

横からの図

高さが激しい図

高さが激しい場合に若干本当にこのように見えるか疑問が残るが、シンプルな波を起こした場合の屈折と反射をray traceで表現してみた。横から覗いた場合、水面だけを取り出して見ているので違和感が残るのかも。

行波について

海面のような進行波については、以下の記事を見つけた。
http://www.rcnp.osaka-u.ac.jp/~hosaka/lecture/gendai_butsuri/waves.pdf
高校2年生で習う物理の内容であることが書かれていて、たしかそんな気もしたなぁと朧気に思い出した。簡単には、進行波


y = f(x-vt)
として表される式で、地点の変数がx-vtの形で表されるようだ。時間ごとに平行移動されていくということですね。上図の式では

y = -1.5 + (sin(x - t) + sin(z - t)) / 8
としています。

反射率について

高さが激しい場合に、水面の高さがmaxのときに屈折した底面が見えたのでこれに違和感を持ち、反射率を見直してみた。結論反射率は関係なかったが、図のような限りある平面の上に起こされた波を実際にみることはまぁないので、この通りになるかどうかわからないですね。CGは現実をトレースするわけではなくきれいに見せれば良いという言い訳をしつつ、理由に気づいたら更新すると思います。
さて、反射率はP偏光とS偏光という2種類の光(電磁波)で計算されるようです。
偏光 - Wikipedia
それぞれの偏光ごとの反射率の計算式は以下などに書かれています。
https://www.chiba-c.ed.jp/funako/fttp_kousin/ssh/reserch/2019/2019_03p3.pdf
9章:反射と屈折の法則
wikiの方のリンクに書かれてますが、これらを50%の配分で足し合わせれば良さそうです。

この式のとおりに計算しましたが、これは反射光と屈折光の配分の話なので、波が高い場合などは関係ないですね。まぁ光を当てたりする場合に綺麗に表現されそうなので良しとしましょう。

水の表現の挑戦

グラフィクスといえば水の表現!!(個人的意見)
ということで vulkan ray tracingで挑戦してみました。
結果はこちら

前回に引き続いて、反射/屈折の部分はガラスキューブと同じです。


ブログということなので、実装するまでで気になったことや工夫した点を主に残そうと思います。

水面の表現

水面は大量のポリゴンで表現しています。図は4.0 * 4.0四方の平面と水面を書いているだけなのですが、平面は三角形のポリゴン2つだけなのに対して、水面は一辺が0.02の三角形を80000ポリゴンで書いています。

波の位置の計算

最初はcompute pipelineで計算していたのですが、おそらく同期とかの部分の実装が甘く、ものすごく遅かったです。上記だと30fps程度ですが、compute pipeline使用時だと1fpsも出てなかったです。上記例では素直にCPUで計算しています。

波の計算式


y = -1.0 + sin(2\pi x) * sin(\pi t) + cos(2\pi z) * sin(\pi t)
としています。vulkanではxz平面に縦軸のyが下向きにあるので、このような式にしています。x方向とz方向で異なる運動の成分を持つようにすれば波っぽくなる気がします。

終わりに

以上、水面の反射/屈折の最も簡単な例を書いてみました。ここから発展させていって、もっときれいな水を表現したいものですね。

ガラスキューブの描画 反射と屈折

Ray tracingの環境が整ったのならば、一度はやってみたいガラスや水の透過、反射、屈折。
ということでvulkan ray tracingでガラスキューブで反射と屈折の描画をしてみたので、その間に調べたことや難しかった点などをまとめておきます。
(一番躓いた点は法線の計算がバグってたという、全く関係のない部分でした;;)
ソースコードの共有というよりはその周りの話がメインとなるので、Qiitaではなくこちらに投稿します。
結果はこんな感じ


f:id:Aqoole_Hateena:20210605201025p:plain
f:id:Aqoole_Hateena:20210605201101p:plain
f:id:Aqoole_Hateena:20210605201123p:plain
f:id:Aqoole_Hateena:20210605201150p:plain

オブジェクトは平面とガラスキューブだけという、非常にシンプルでray tracing学習者丸出しの内容です。

さて、このトレースではガラスキューブの入射光から、反射光と屈折光という2種類の計算を行っているのですが、何が難しいかと言うと
・シェーダ内で関数traceRayEXT()は再帰呼び出しできない(現状AMD GPUのみ?)
・GLSLでは一般の関数でも再帰呼び出しはできない(or 非推奨)
・しかしながら、ガラスに当たる光の計算は1種類の入射光から反射光と屈折光の2種類が生まれる

はい、そういうことです。
ガラス面に当たるたびに、ひとつの入射光から2つの光が生まれるので、2の指数の速さで光の本数が増えていき、コードの実装が大変です。例えばC++だとこのようなケースでは再帰関数で書いておいて、進んだり戻ったりを繰り返しトレースできそうなものですが、GLSLでかくシェーダになると再帰関数は使えないそうです。(ソースは若干古いです)

https://stackoverflow.com/questions/43601521/recursion-in-glsl-prohibited/43601658

GLSLで再帰関数が使えないのならtraceRayEXT関数を再帰的に使えばいいじゃない、という発想になります。つまりtraceRayEXT()でヒットした場合の呼び出し先のシェーダ内でさらに同じtraceRayEXT()を使うということですね。しかし結論的に、これはどうもできないようです。微妙な言い方になるのは、下記のNVIDIAの記事だと「遅くはなるが動く」と言っているようですが、私のAMDGPUだとVkPhysicalDeviceRayTracingPipelinePropertiesKHR内のmaxRayRecursionDepthという値を調べると1だったので、1回しか呼べないということになっています。つまり、再帰呼び出しできません。

https://github.com/nvpro-samples/vk_raytracing_tutorial_KHR/tree/master/ray_tracing_reflections

仕方がないので、再帰関数でスマートに書くことは諦めてトレースするレイの本数を直打ちで決めてトレースしました。最終的にはトレースする光線の数は2の指数の和になります。
この例ではdepth = 6, つまり6回光がガラス面にあたって反射/屈折した場合のトレースをしているので、1ピクセルあたり


\begin{eqnarray}
\sum_{k=1}^{6}2^k &=& 2^7 - 2 \\
&=& 126\
\end{eqnarray}

の光をトレースしています。
あらかじめ

vec3[126] colors;

に126本の光線のトレースをしてヒット先の色を格納した後、反射光の強度や全反射の有無にあわせて

/*
for(uint i = 62; i < 126; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}
*/
for(uint i = 30; i < 62; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}

for(uint i = 14; i < 30; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}

for(uint i = 6; i < 14; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}
for(uint i = 2; i < 6; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}
for(uint i = 0; i < 2; i++)
{
  uint reflectIndex = 2 * i + 2;
  uint refractIndex = 2 * i + 3;
  colors[i] = ColorBlend(reflection, colors[i], isPlane[i], isMiss[reflectIndex], colors[reflectIndex], isMiss[refractIndex], colors[refractIndex], isAllReflect[refractIndex]);
}
return ColorBlend(reflection, surfaceColor, false, isMiss[0], colors[0], isMiss[1], colors[1], isAllReflect[1]);
vec3 ColorBlend(float reflection, vec3 ownColor, bool isPlane, bool reflectMiss, vec3 reflectColor, bool refractMiss, vec3 refractColor, bool isAllReflect)
{
  if(isPlane)
    return ownColor;
  if(reflectMiss && refractMiss)
    return vec3(0.0, 0.0, 1.0);
  if(reflectMiss)
  {
    if(isAllReflect)
      return ownColor;
    else 
      return mix(refractColor, ownColor, 0.1);
  }
  if(refractMiss)
    return mix(reflectColor, ownColor, 0.1);
  //both hit
  if(isAllReflect)
    return mix(reflectColor, ownColor, 0.1);
  else
  {
    vec3 traceColor = mix(refractColor,reflectColor,reflection);
    return mix(traceColor,ownColor,0.1);
  }
}

などとして色をブレンドしていきます。
この場合だとガラスから抜け出て平面にヒットした場合でも構わずトレースし続けているので改善の余地はありますが、予めトレースする光線の数を決めておくことに変わりはありません。GLSLではmalloc/newのように動的にメモリを確保する手段がないはずなので。
最初の/* */でコメントアウトしている部分はその次のdepthまでトレースした場合の処理ですが、ここまで行うと、どうやらvulkanの側でpipelineを作成する際にVK_ERROR_UNKNOWN = -13でエラーしていましたので、その一つ前のdepthで諦めています。fps的にはdepthが6の時点で400〜700fpsなので、もう少し余裕はあったかなと思ってます。

colors[126]にヒットした先の色が取ってこれたとして、次のColorBlend()の処理で色を合成しています。この処理はdepthが深いところの色を決めてから、その色でもって浅いところの色が決まっていきます。
ルールとしては
・平面にヒットした場合は、前の結果(よりdepthが大きい位置)にかかわらず平面の色を返す
・反射/屈折がミス(どこにもヒットしない)した場合、どちらかの色をガラス面とミックスして返す
・全反射が起きた場合、屈折光の色をガラス面とミックスして返す
・反射/屈折の両方がヒットしている場合は、反射率だけ反射光と屈折光の色を合成し、それをガラス面とミックスして返す
ガラス面とミックスする部分は、色の綺麗さなどに合わせてお好みかもしれません。
反射率は以下のサイトを参考に計算しています。
http://skomo.o.oo7.jp/f17/hp17_9.htm
https://ja.wikipedia.org/wiki/%E5%8F%8D%E5%B0%84%E7%8E%87

イデアとしてはたったこれだけなんですが、試行錯誤だったり何だったりで非常に時間がかかりました。とりあえず結果が得られたということで、満足しています。