※このブログはVladimir Agafonkinが執筆した記事の翻訳です。

以下、私＝ウラジミール氏

WebGLを使った風のシミュレーションのデモをご覧ください。それでは、実際にどのように動いているのかを見てみましょう。

‍

お恥ずかしい話ですが、Mapbox社で働いていたこの数年間、私はOpenGL/WebGLを使った直接的なプログラミングを避けていました。その理由は、OpenGLのAPIや用語があまりにも複雑で、わかりにくい印象があったので、思い切って飛び込んでみることができませんでした。ミップマップ、ブレンド関数などの用語を聞いただけで、不安な気持ちになったものです。

私はついに自分の恐怖に立ち向かい、WebGLを使って自力で何か非日常的なものを作ってみることにしました。2Dの風シミュレーションはとても良い機会でした。便利で、視覚的に美しく、挑戦的でありながら、まだ達成可能な範囲だと感じたのです。そして、WebGLは、実際見た目よりもずっと怖くないことに驚きました。

CPUによる風の可視化

風のビジュアライゼーションの例はネット上にたくさんありますが、最も人気があり影響力のあるものは、キャメロン・ベッカリオ氏の有名なプロジェクトであるearth.nullschool.netです。このプロジェクト自体はオープンソースではありませんが、古いオープンソースのバージョンがあり、他のほとんどの実装がそのコードをベースにしています。注目すべきオープンソースの派生として、Esri Wind-JSがあります。この技術を使った人気の高い気象サービスには、WindyやVentuSkyなどがあります。

‍

一般的に、ブラウザでのこのような可視化は、Canvas 2D APIに依存しており、大まかには次のように動作します。

1.画面上にランダムな粒子位置の配列を生成し、粒子を描画する。

2.各粒子について、風のデータを照会して現在の位置での粒子速度を取得し、それに応じて移動させる。

3.粒子の一部をランダムな位置にリセットする。これにより、風が吹き飛ばされた場所が完全に空にならないようにする。

4.現在の画面を少しフェードして、新しく配置された粒子を上に描画します。

‍

これにはパフォーマンス上の制限があります。

• 風の粒子数は少なくする必要がある（例：earth.nullschool.netは5k程度）。
• データ処理は高価でCPU側で行われるため、データやビューを更新するたびに大きな遅延が発生する（例：地球では約2秒）。

さらに、MapboxのようなWebGLベースのインタラクティブな地図の一部として統合するには、フレームごとにCanvas要素のピクセルコンテンツをGPUにアップロードする必要があり、パフォーマンスが大幅に低下してしまいます。

‍

私は、WebGLを使ってGPU側で完全なロジックを再実装する方法を探していました。そうすれば、高速で、何百万もの粒子を描くことができ、大きなパフォーマンスの低下なしにMapbox GLマップに統合することが可能になります。幸運なことに、クリス・ ウェロンス氏によるWebGLでの粒子物理学に関する素晴らしいチュートリアルを偶然見つけ、風のビジュアライゼーションにも同じアプローチが使えることに気付きました。

OpenGLの基本

APIや用語が混乱しているため、OpenGLグラフィックスプログラミングを学ぶのは非常に難しいのですが、表面上は非常にシンプルなコンセプトです。実用的な定義を紹介しましょう。

OpenGLは、三角形を効率的に描画するための2D APIを提供します。

つまり、GLで行うことは、基本的に三角形を描くことだけなのです。難しいという印象は、これを行うために必要な様々な数学やアルゴリズムから来ています。点や基本的な線（スムージングや丸いつなぎ目/キャップなし）も描けますが、これらはほとんど使われません。

OpenGLでは、GPUが直接実行するプログラムを記述するために、C言語に似た特別な言語であるGLSLが用意されています。各プログラムはシェーダーと呼ばれる2つの部分（バーテックスシェーダーとフラグメントシェーダー）に分かれています。

バーテックスシェーダは、座標を変換するためのコードを提供します。例えば、三角形の座標を2倍にして、三角形が2倍の大きさに見えるようにするなどです。描画時にOpenGLに渡す座標ごとに1回実行されます。基本的な例です。

attribute vec2 coord;
void main() {
   gl_Position = vec4(2.0 * coord, 0, 1);
}

フラグメントシェーダーは、描画された各ピクセルの色を決定するためのコードを提供します。いろいろな計算をすることができますが、最終的には「三角形の現在のピクセルを緑で描画する」というようなことになります。

void main() {
   gl_FragColor = vec4(0, 1, 0, 1);
}

バーテックスシェーダーでもフラグメントシェーダーでも、パラメータとして画像（テクスチャと呼びます）を追加し、その画像の任意の場所のピクセルカラーを調べることができるという優れた機能があります。風のビジュアライゼーションでは、この機能を多用します。

フラグメントシェーダのコード実行は超並列かつハードウェアアクセラレーションを多用しているため、通常はCPUでの同等の計算よりも飛躍的に高速です。

風のデータの取得

米国国立気象局は、6時間ごとに地球全体の気象データ（GFS）を、緯度・経度のグリッドとそれに関連する値（風速など）の形式で発表しています。このデータはGRIBと呼ばれる特殊なバイナリ形式でエンコードされており、特殊なツールを使って人間が読めるJSONに解析することができます。

私は、風のデータをダウンロードして、風速をRGBカラーでエンコードしたシンプルなPNG画像に変換するいくつかの小さなスクリプトを書きました（各ピクセルの水平方向の速度は赤、垂直方向の速度は緑）。これは次のようなものです。

‍

もっと高解像度のものもありますが（2倍、4倍）、360×180グリッドであれば、低ズームでの視覚化には十分です。この種のデータにはPNG圧縮が非常に適しており、上の画像の重さは通常約80KBです。

GPUによる粒子の移動

既存の風のビジュアライゼーションでは、粒子の状態をJavaScriptの配列に格納していました。では、この状態をGPU側でどのように保存・操作すればよいのでしょうか。OpenGL ES 3.1およびそれに相当するWebGL 2.0の仕様では、コンピュートシェーダーと呼ばれるGLの新機能により、任意のデータに対してシェーダコードを実行することができます（レンダリングは行いません）。しかし、残念ながら、ブラウザやモバイル機器などの新仕様への対応はわずかで、現実的な選択肢は「テクスチャ」しかありませんでした。

OpenGLでは、スクリーンだけでなく、フレームバッファと呼ばれる概念を使って、テクスチャにも描画することができます。そのため、粒子の位置を画像のRGBAカラーとしてエンコードし、それをGPUにロードし、フラグメントシェーダで風速に基づいて新しい位置を計算し、再びRGBAカラーにエンコードして、新しい画像に描画することができます。

XとYの両方に十分な精度を持たせるために、それぞれの成分をRGとBAの2バイトに格納し、65536個の異なる値の範囲を確保しています。

500×500の画像に25万個の粒子を配置し、フラグメントシェーダで全ての粒子を動かしてみます。その結果、画像は次のようになります。

‍

フラグメントシェーダーでは、RGBAから位置をデコードしたり、エンコードしたりしています。

// lookup particle pixel color
vec4 color = texture2D(u_particles, v_tex_pos);
// decode particle position (x, y) from pixel RGBA color
vec2 pos = vec2(
   color.r / 255.0 + color.b,
   color.g / 255.0 + color.a);
... // move the position
// encode the position back into RGBA
gl_FragColor = vec4(
   fract(pos * 255.0),
   floor(pos * 255.0) / 255.0);

次のフレームでは、この新しい画像を現在の状態とし、新しい状態をもう一方の画像に描画するなど、フレームごとに2つの画像を入れ替えていきます。このように、2つの粒子状態テクスチャを使用することで、風のシミュレーションロジックをすべてGPUに移すことができます。

この方法は非常に高速で、ブラウザ上で毎秒60回、5,000個の粒子を更新するだけだったのが、いきなり100万個の粒子を処理できるようになります。

留意点としては、極付近では赤道上の粒子に比べてX軸方向の移動速度が速くなることが挙げられますが、これは経度の数が同じであれば距離がずっと短くなるためです。これを考慮して、以下のシェーダーコードを作成しました。

float distortion = cos(radians(pos.y * 180.0 - 90.0));
// move the particle by (velocity.x / distortion, velocity.y)

粒子の描画

先に述べたように、三角形に加えて基本的な点も描くことができます。あまり使われることはありませんが、今回のような1ピクセルの粒子には最適です。

各粒子を描画するには、バーテックスシェーダーで粒子の状態を表すテクスチャのピクセルカラーを調べて位置を決定し、フラグメントシェーダーで風のテクスチャから現在の速度を調べて粒子の色を決定し、最後に綺麗なカラーグラデーションにマッピングします（信頼できるColorBrewer2から色を選びました）。この時点では、以下のようになります。

すこし空白が目立ちますが、いい感じですね。しかし、粒子の動きだけでは、風向きを感じ取ることができません。粒子の軌跡を追加する必要があります。

粒子の軌跡を描く

最初に試した軌跡の描き方は、画面の状態をフレーム間で持続させるWebGLのオプション「preserveDrawingBuffer」を使うことで、粒子の動きに合わせてフレームごとに何度も重ねて描くことができます。このWebGLの機能は高いパフォーマンスを示すため、多くのWebGLの記事では使用を推奨しています。

その代わり、粒子の状態を表すテクスチャと同じように、粒子をテクスチャに描画し（テクスチャはスクリーンに描画されます）、そのテクスチャを次のフレームの背景として使用し（少し薄暗くします）、入力/ターゲットのテクスチャをフレームごとに入れ替えます。この方法の利点は、パフォーマンスの向上だけでなく、（preserveDrawingBufferに相当するものがない）ネイティブコードに直接移植することができることです。

補間された風のルックアップ

‍

‍

風のデータは、緯度・経度グリッド上の特定のポイント、例えば(50,30),(51,30),(50,31),(51,31)の地理的ポイントの値を持っています。では、任意の中間値、たとえば(50.123,30.744)を得るにはどうすればよいのでしょうか。

OpenGLでは、テクスチャの色を調べるときに、自由に補間を行うことができます。しかし、この方法では、ブロック状のピクセルパターンになってしまいます。以下は、風のテクスチャを拡大したときのアーチファクトの例です。

‍

幸いなことに、各ウィンドプローブ内の4つの隣接ピクセルを調べ、フラグメントシェーダー内のネイティブな補間処理に加えて、手動でバイリニア補間処理を行うことで、アーチファクトを滑らかにすることができます。修正に多少の時間はかかりますが、アーチファクトが修正され、より滑らかな風の可視化が可能になります。このテクニックを使った同じエリアの写真です。

GPU上の擬似乱数ジェネレータ

GPUに実装するには、粒子の位置をランダムにリセットするという、トリッキーなロジックを行う必要があります。これがないと、膨大な数の風の粒子があっても、時間が経つと風が吹いた場所が空になってしまうため、画面上の数行にしかなりません。

問題は、シェーダーには乱数発生器がないことです。粒子をリセットする必要があるかどうかを、どうやってランダムに判断するのでしょうか？

私はStackOverflowで解決策を見つけました。数字のペアを入力として受け取る、疑似乱数生成用のGLSL関数です。

float rand(const vec2 co) {.
   float t = dot(vec2(12.9898, 78.233), co);
   return fract(sin(t) * (4375.85453 + t))。
}

この不思議な関数は、sinの結果が大きな値で大きく変化することに依存しています。そこで、次のようなことができます。

if (rand(some_numbers) > 0.99)
   reset_particle_position();

問題は、生成された値が画面全体で均一になり、奇妙なパターンを示さないように、十分に「ランダム」であるような各粒子の入力を選ぶことです。

現在の粒子の位置をシードとして使用することは、同じ粒子の位置では常に同じ乱数が生成されるため、同じエリアで消える粒子が出てくるなど、完璧ではありません。

また、粒子の位置を状態テクスチャに使用しても、同じ粒子が常に消えてしまうため、うまくいきません。

結局、粒子の位置と状態の位置の両方に依存して、フレームごとに計算されたランダムな値をシェーダーに渡すことにしました。

vec2 seed = (pos + v_tex_pos) * u_rand_seed;

しかし、もう一つの問題があります。粒子の速度が非常に速い場所は、風があまり吹いていない場所に比べてはるかに密集して見えます。粒子が速い場合は粒子のリセットレートを上げることで、少しバランスを取ることができます。

float dropRate = u_drop_rate + speed_t * u_drop_rate_bump;

ここでの speed_t は相対的な速度の値 (0 から 1) で、u_drop_rate と u_drop_rate_bump は最終的なビジュアライゼーションで微調整できるパラメータです。結果にどのような影響を与えるかの例を示します。

‍

今後の展開

結果として、100万個の粒子を60fpsでレンダリングできる、完全にGPUを搭載した風のビジュアライゼーションが完成しました。デモのスライダーで遊んでみたり、最終的なコードをチェックしてみてください。

次のステップは、これをライブマップに統合して探索できるようにすることです。少しずつ進めていますが、ライブデモを公開するほどではありません。画質が悪くて申し訳ありませんが、少しだけ見てみましょう。

お読みいただきありがとうございました。

グラフィックスプログラミングに関する私のくだらない質問に辛抱強く答えてくれて、たくさんの貴重なヒントを与えてくれたMapboxチームメイトには、心より感謝しています。