【WebGPU】Compute Shader で Curl Noise を計算してパーティクルを動かす

🎄この記事は【カヤック】面白法人グループ Advent Calendar 2024の18日目の記事です 🎄

こんにちは！ハイパーカジュアルゲームチーム・エンジニアの深澤です。

WebGPU の Compute Shader で Curl Noise を計算し、パーティクルの位置を更新してみました。
スクショは、MacBookPro M1で100万個のパーティクルを動かしたものです。

画像をクリックするとデモに飛びます。
WebGPU で実装しているため、Chromeのみの動作となります。

デモURL: https://takumifukasawa.github.io/webgpu-particle-compute-shader-curl-noise-demo/

リポジトリ github.com

動機

たとえばこの動画のように、大量のパーティクルが一つ一つまるで風に揺れるように・流体のように動く様子はかっこいいですよね！

www.youtube.com

このような動きを、流体の方程式に沿ってちゃんとシミュレーションをするとどうしても大変な計算負荷になります。
リアルタイムレンダリングとなると、初期位置や速度の決定から、動いている途中の位置・速度を毎フレーム更新することが必要になります。
数十、数百のパーティクルぐらいなら耐えられるかもしれませんが、数千、数万と数が多くなるにしたがってより実現が難しくなっていくはずです。

そこで、Curl Noise という手法を使うことにより流体のシミュレーションをせずとも流体に近い動きが実現できます。

自分で Curl Noise の実装をしたことがなかったので、アドベントカレンダーを機に試してみることにしました。

せっかくなので、Chrome で標準で使用できるようになった WebGPU と Compute Shader を使っていきます。
WebGL でも TransformFeedback を使うことで同等のことは可能です。

Curl Noise

※理解が間違っている部分などありましたらすいません！

概要

SIGGRAPH 2007 で発表された、疑似流体アルゴリズムの発表の中に登場します。論文/資料はこちらになります。

https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph2007-curlnoise.pdf

この手法が開発された理由として、既存の流体運動の方程式によるシュミレーションはリアルな挙動を表現するが、計算速度の問題やパラメーターによる制御のコントロールが難しいことが発端にあったようです。
ベクトルのポテンシャル場が2次元ないし3次元空間に構築されており、速度をベクトルのポテンシャル場から求める、という方法です。

まず、ベクトルのポテンシャル場とは空間にベクトルが散らばっているイメージです。以下、ベクトル場と呼びます。
風を例にすると、屋外では風が横に吹いて樹木が揺れたり、木の葉が上に舞い上がったりしますよね。つまり、とある空間の中に様々な方向の力があちこちに存在している、というイメージです。

2次元の場合のイメージ。平面上にベクトルが散らばっている

- ベクトル場の中から、とある地点のベクトルの勾配（後述）を取り出す
- そのベクトルを回転させた値（ベクトル解析における回転）を速度として扱う
- 現在の位置を更新していくことで渦のような動きを得ることができる
というのが Curl Noise の核のようです。

特徴

Curl Noise の特徴として「完全に非圧縮性」という性質があります。非圧縮性というのは、密度が一定で、発散が0になるということのようです。速度場においては、入力と出力の速度が一定、という理解です。

計算

∇（= ナブラ）はベクトル微分演算子のことで、この文脈の中ではとても大雑把にいうと「勾配」として扱います。
計算には、現在の位置からx,y,zを少しずらした位置のベクトルを取得し、勾配をとったものです。

そして、ナブラの外積をとってベクトルを回転させたものが Curl Noise の計算結果となります。
各参考記事によるとこういう計算をすることになるようです。
（すいませんこのあたりは知識が特になく、理解が甘いです）

また、回転したベクトルと元のベクトルの内積をとったものが「発散」です。
非圧縮性の性質から発散が0になる = 内積が0 = 90度なので「ベクトル解析において90度回転させたものが Curl Noise の結果」と言い換えることもできそうです。

...つまり細かい計算は分かっていない部分が多々あるのですが、
ベクトル場を生成して、勾配をとり、ベクトルを回転させ、それを速度として捉える、ということなのかなと思っています。

本記事で触れること / 触れないこと

論文では境界制限（衝突,範囲制限）や粘性についても触れていますが、今回は単純に速度を求めることのみに着目します。

全体の流れ

速度の計算を Compute Shader で行うので座標の更新も Compute Shader で行っていきます。

大量にパーティクルを出すために、板ポリを GPU Instancing させます。
板ポリにする理由は、Point の描画ではなく、大きさを持ったビルボードとして常にカメラを向かせるような描画をさせたかったためです。

パーティクルのシェーダー内でのインスタンシング情報の読み取りですが、Compute Shader で GPU に紐づいたBufferを更新する方針なので、そのBuffer を共有・参照すれば計算をGPUで完結することができるので、CPU への読み書きや Buffer のバインドし直しなどが発生せず効率が良くなりそうです。

パーティクルを描画するためのRenderPipelineは以下のような構造になります。

// positionとvelocityをvec4<f32>な値として詰め込む
const instanceDataArray = new Float32Array([
    0, 0, 0, 0, // インスタンス=0のposition
    0, 0, 0, 0, // インスタンス=0のvelocity
    ...
]);

const particleInstancesBuffer = gDevice.createBuffer({
    size: instanceDataArray.byteLength,
    usage:
        GPUBufferUsage.STORAGE |
        GPUBufferUsage.VERTEX |
        GPUBufferUsage.COPY_DST,
    mappedAtCreation: true
});
new Float32Array(particleInstancesBuffer.getMappedRange()).set(instanceDataArray);
particleInstancesBuffer.unmap();

...

const particleRenderPipeline = gDevice.createRenderPipeline({
    layout: 'auto',
    vertex: {
        module: gDevice.createShaderModule({
            code: vertWGSL, // 頂点シェーダー
        }),
        entryPoint: 'main',
        buffers: [
            // buffer = 0
            {
                arrayStride: quadVertexSize, // UVなど板ポリの情報で使う頂点データの大きさ
                stepMode: 'vertex',
                attributes: [
                    // 頂点が持つデータ。uvなど任意に設定
                    {
                        shaderLocation: 0, // @location(0)
                        offset: quadPositionOffset,
                        format: 'float32x4'
                    },
                    {
                        shaderLocation: 1, // @location(1)
                        offset: quadColorOffset,
                        format: 'float32x4'
                    },
                    {
                        shaderLocation: 2, // @location(2)
                        offset: quadUVOffset,
                        format: 'float32x2'
                    }
                ]
            },
            // buffer = 1
            {
                arrayStride: instancingVertexSize, // インスタンシングごとの頂点データの大きさ
                stepMode: 'instance', // Bufferをインスタンシングとして扱う
                attributes: [
                    // インスタンシングごとのデータ。今回はpositionとvelocityを格納したBufferを扱う
                    {
                        shaderLocation: 3, // @location(3)
                        offset: 0,
                        format: 'float32x4'
                    },
                    {
                        shaderLocation: 4, // @location(4)
                        offset: 4 * 4,
                        format: 'float32x4'
                    }
                ]
            }
        ]
    },
    ...
});

頂点シェーダーはこのような構造になっています。

struct Particle {
    position: vec4<f32>,
    velocity: vec4<f32>
};

@group(1) @binding(0) 
var<storage, read> instanceData : array<Particle>;
  
...  
@vertex 
fn main(
    // buffer=0の頂点データを展開
    @location(0) position : vec4<f32>, // ビルボードのオフセット方向
    @location(1) color : vec4<f32>,
    // buffer=1のインスタンスデータを展開
    @location(2) uv: vec2<f32>,
    @location(3) instancePosition: vec4<f32>,
    @location(4) instanceVelocity: vec4<f32>
) -> VertexOutput {
     ...

仮に、板ポリの描画用の情報が入ったBufferを particleVerticesBuffer、インスタンシングのデータが入った Buffer を particleInstancesBuffer とすると、以下のようにBufferを登録することができます。

...

const renderPassEncoder = renderCommandEncoder.beginRenderPass(renderPassDescriptor);
renderPassEncoder.setPipeline(particleRenderPipeline);

... uniforms の設定など

renderPassEncoder.setVertexBuffer(0, particleVerticesBuffer); // UVなど板ポリの情報が入った頂点のBufferをindex=0に登録
renderPassEncoder.setVertexBuffer(1, particleInstancesBuffer); // 全てのインスタンシングのデータが詰まったBufferをindex=1に登録

// インデックスを使ったインスタンシング描画
renderPassEncoder.setIndexBuffer(quadIndicesBuffer, 'uint16');
renderPassEncoder.drawIndexed(quadIndexArray.length, instanceNum);

そして、この particleInstancesBuffer を Compute Shader で直接更新していきます。

Compute Shader で実装

まずは Compute Shader の実行関数から見ていきましょう。

curlNoiseという関数の結果から速度を算出しています。

struct Instance {
    position: vec4<f32>,
    velocity: vec4<f32>
};

@group(0) @binding(0)
var<uniform> uniforms : Uniforms;
            
@group(0) @binding(1)
var<storage, read_write> input : array<Instance>; // インスタンスのバッファを読み込み、書き込み先にも使う

...

// jsからシェーダーを template literal として書き、workGroupSizeを埋め込む
// 引数が一つの場合、一次元になる
@compute @workgroup_size(${workgroupSize}) 

fn main(
    @builtin(global_invocation_id) global_invocation_id: vec3<u32>
) {
    let id = global_invocation_id.x; // 一次元なので、global_invocation_id をそのままインスタンスIDとして扱うことができる
    if(id < arrayLength(&input)) {
        let instance = input[id];
        let currentPosition = instance.position.xyz;
        let currentVelocity = instance.velocity.xyz;
        let force = curlNoise(currentPosition.xyz * uniforms.noiseScale) - currentVelocity;
        let newVelocity = force * uniforms.speed * uniforms.deltaTime;
        let newPosition = currentPosition + newVelocity;
        input[id].position = vec4<f32>(newPosition.xyz, 1.);
        input[id].velocity = vec4<f32>(newVelocity.xyz, 1.);
    }
}

javascript 側はこのようになっています。

...

const computeBindGroup = gDevice.createBindGroup({
    layout: computePipeline.getBindGroupLayout(0), // @group(0)
    entries: [
        {
            binding: 0, // @binding(0)
            resource: {
                buffer: computeUniformBuffer
            },
        },
        {
            binding: 1, // @binding(1)
            resource: {
                buffer: particleInstancesBuffer
            },
        },
    ]
});

...

const computeCommandEncoder = gDevice.createCommandEncoder();
const computePassEncoder = computeCommandEncoder.beginComputePass();
computePassEncoder.setPipeline(computePipeline);
computePassEncoder.setBindGroup(0, computeBindGroup);
computePassEncoder.dispatchWorkgroups(dispatchGroupSize); // 引数が一個なら一次元の実行
computePassEncoder.end();
gDevice.queue.submit([computeCommandEncoder.finish()]);

今回の実装では、workGroupとdispatchがどちらも一次元なので、
実行されるスレッド数は workgroup_size * dispatchWorkgroupsの数 となります。
この数が描画したいパーティクルの数と等しくなるように扱うことで、表示されているパーティクルをすべて更新することが可能となります。

それぞれ上限が決まっており、論理デバイスから値を取得することができます。
workGroupSizeは3次元を、dispatchWorkgroup は1つの数値を取得することができます（なので、xyzすべて同じ値だと予想します）。

以下、参考コードです。

const gAdapter = await navigator.gpu.requestAdapter(); // 物理デバイス
const gDevice = await gAdapter.requestDevice(); // 論理デバイス

// デバイスの制限値を取得
const limits = gDevice.limits;

// 手元のMacbookPro M1では256,64,64
console.log("Max Compute Workgroup Size:", {
    x: limits.maxComputeWorkgroupSizeX,
    y: limits.maxComputeWorkgroupSizeY,
    z: limits.maxComputeWorkgroupSizeZ,
});

// 手元のMacbookPro M1では65536
console.log("Max Compute Workgroups Per Dimension:",
    limits.maxComputeWorkgroupsPerDimension,
);

手元のMacbookProではそれぞれ x:256,y:64,z:64、65536でした。
つまり、workGroupSizeとdispatchWorkgroupを一次元で扱おうとする場合は 256 * 65536 = 16777216 のスレッドを同時に取り扱うことができる、ということになります。

Curl Noise の計算

論文ではベクトル場の生成に Perlin Noise が使われています。

今回の実装では Simplex Noise を使うこととします。
Simple Noise の実装はこちらを元にさせていただきます。

fn curlNoise(position: vec3<f32>) -> vec3<f32> {
    let eps = 0.001;
    let eps2 = 2.0 * eps;
    let invEps2 = 1.0 / eps2;
    let dx = vec3<f32>(eps, 0.0, 0.0);
    let dy = vec3<f32>(0.0, eps, 0.0);
    let dz = vec3<f32>(0.0, 0.0, eps);
     
    // 勾配算出のため、epsだけずらした地点のベクトルを算出  
    let px0 = snoise3(position - dx);
    let px1 = snoise3(position + dx);
    let py0 = snoise3(position - dy);
    let py1 = snoise3(position + dy);
    let pz0 = snoise3(position - dz);
    let pz1 = snoise3(position + dz);

    // 回転
    let x = (py1.z - py0.z) - (pz1.y - pz0.y);
    let y = (pz1.x - pz0.x) - (px1.z - px0.z);
    let z = (px1.y - px0.y) - (py1.x - py0.x);

    return vec3<f32>(x, y, z) * invEps2;
}

入力に使う座標からxyzを少しずらした位置のベクトルをベクトル場から取り出し、勾配を求め、ベクトルの回転をかける計算です。

simplex noise の調整

simplex noise の元の実装が vec3 の引数を元に float を返す実装になっているので、x,y,zそれぞれにsnoiseの結果を割り当てることにします。

このときポイントとなるのが、非連続になりづらいように「大きくオフセットした」ノイズの値にすることのようです。
おそらく、元のノイズのパターンの影響を受けづらいようにするためなのかなと思います。

fn snoise3(v: vec3<f32>) -> vec3<f32> {
    return vec3<f32>(
        snoise(v),
        snoise(v + vec3<f32>(100., 200., 300.)), // 任意のオフセット量を適用
        snoise(v + vec3<f32>(400., 500., 600.))
    );
}