ローレンツアトラクタを描画してみる

UnityでGPGPUシリーズはその1~その4まで数値計算と描画を切り分けて、数値計算に特化した内容で記事をまとめてきた。ただやはり可視化してみないと何を計算しているのかわからないことも多いので、ここでは最低限の可視化について触れようと思う。
タイトルにあるようDrawProceduralNowを使えば結構楽にGPUインスタンシング(であってる?)ができるのでまとめてみた。

可視化の題材はローレンツアトラクタ
詳しくはググってもらうとして、これはある計算をすることで見た目が3Dのこんな感じの絵がでるようになる。

lorenz

適当なx,y,z座標の初期値を与えると、あとは規則に沿ってこういう軌道を描き、なんかちょっと面白い絵ができる。

これは普通CPUで計算するような題材だが、計算も簡単なのでGPU上で座標計算を行い、可視化までGPU内で完結させてしまおうと思う。ComputeBufferにx,y,z座標を格納し、その座標に1ピクセルのドットを描画する、というのをまずは目指す。

実装

今回必要なコードはたった3つ
・C#コード
・ComputeShaderコード(座標計算用)
・SurfaceShader.shader(可視化用コード)

でこれに加えて
・空のGameObject


まず空のGameObjectを作成し
Assetsフォルダ内にScriptsフォルダを作成
・右クリック→Create→C# script→「Main」
Assetsフォルダ内にShaderフォルダを作成
・右クリック→Create→Shader→Compute Shader→「ComputeShader」
・右クリック→Create→Shader→Standard Surface Shader→「SurfaceShader」

以下のソースをコピーして
unitygamen
こんな感じにできれば完成。


Main.cs
using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class Main : MonoBehaviour
{
    //計算するコンピュートシェーダー
    [SerializeField]
    ComputeShader cshader;

    // ドットをレンダリングするシェーダー
    [SerializeField]
    Shader renderingShader;
    Material renderingShader_Material;

    ComputeBuffer posBuffer;

    //[SerializeField]
    int particleNum = 1024 * 1024;

    int kernel_Lorenz;
    int kernel_BufInit;

    void Start()
    {
        renderingShader_Material = new Material(renderingShader);

        posBuffer = new ComputeBuffer(particleNum, sizeof(double) * 3);//x,y,z成分で3つ
        kernel_Lorenz = cshader.FindKernel("Lorenz");
        kernel_BufInit = cshader.FindKernel("BufInit");
        cshader.SetBuffer(kernel_Lorenz, "posBuffer", posBuffer);
        cshader.SetBuffer(kernel_BufInit, "posBuffer", posBuffer);
        
        //初期化もGPUで
        cshader.Dispatch(kernel_BufInit, particleNum / 64, 1, 1);

        // GPUバッファをマテリアルに設定
        renderingShader_Material.SetBuffer("posBuffer", posBuffer);
    }

    // Update is called once per frame
    void Update()
    {
        cshader.Dispatch(kernel_Lorenz, particleNum / 64, 1, 1);
    }

    void OnRenderObject()
    {
        // レンダリングを開始
        renderingShader_Material.SetPass(0);
        // n個のオブジェクトをレンダリング
        Graphics.DrawProceduralNow(MeshTopology.Points, particleNum);
    }

    private void OnDestroy()
    {
        //解放
        posBuffer.Release();
    }
}


ComputeShader.compute
#pragma kernel Lorenz
#pragma kernel BufInit

//OpenCLに移植するときはfloat3,double3が使えないので注意
RWStructuredBuffer<double3> posBuffer;
#define DELTA_TIME (0.0001)
#define LOOPNUM (128)


//簡易ランダム
uint wang_hash(uint seed)
{
	seed = (seed ^ 61) ^ (seed >> 16);
	seed *= 9;
	seed = seed ^ (seed >> 4);
	seed *= 0x27d4eb2d;
	seed = seed ^ (seed >> 15);
	return seed;
}

//ランダム位置に初期値を置く
[numthreads(64, 1, 1)]
void BufInit(uint id : SV_DispatchThreadID)
{
	double3 pos;
	uint r = wang_hash(id * 1847483629);
	pos.x = 0.005 * (r % 1024) - 1.5;
	r /= 1024;
	pos.y = 0.005 * (r % 1024) - 3.5;
	r /= 1024;
	pos.z = 0.005 * (r % 1024) + 23.0;
	posBuffer[id] = pos;
}

[numthreads(64,1,1)]
void Lorenz(uint id : SV_DispatchThreadID)
{
	double3 pos = posBuffer[id];
	double3 pos_after;
	for (int i = 0; i < LOOPNUM; i++) 
	{
		pos_after.x = pos.x + DELTA_TIME * (-10.0 * pos.x + 10.0 * pos.y);
		pos_after.y = pos.y + DELTA_TIME * (-pos.x * pos.z + 28.0 * pos.x - pos.y);
		pos_after.z = pos.z + DELTA_TIME * (pos.x * pos.y - 8.0 * pos.z / 3.0);
		pos = pos_after;
	}

	posBuffer[id] = pos_after;
}


SurfaceShader.shader
Shader "Custom/SurfaceShader" {
	SubShader{
		ZWrite Off
		Blend One One//加算合成

		Pass {
			CGPROGRAM
			// シェーダーモデルは5.0を指定
			#pragma target 5.0

			#pragma vertex vert
			//#pragma geometry geom
			#pragma fragment frag
			#include "UnityCG.cginc"
	
	StructuredBuffer<double3> posBuffer;

	// 頂点シェーダからの出力
	struct VSOut {
		float4 pos : SV_POSITION;
	};

	// 頂点シェーダ
	VSOut vert(uint id : SV_VertexID)
	{
		// idを元に、ドットの情報を取得
VSOut output; output.pos = mul(UNITY_MATRIX_VP, float4(float3(posBuffer[id]), 1)); return output; } // ピクセルシェーダー fixed4 frag(VSOut i) : COLOR { return float4(0.0154,0.05,0.3,1); } ENDCG } } }

参考
[Unity]コンピュートシェーダ(GPGPU)で1万個のパーティクルを動かす

結果

実行するとこんな感じにきれいな絵が出来上がる。
Webp.net-gifmaker

このgif絵はCameraをわざわざ回転させてるので実際の実行画面とはちょっと違うが。
Camera回転版はGithub


解説

今回の新しい点はSurfaceShader.shaderでComputeBufferを受け取っているというところ。その紐づけはC#側で行う。
renderingShader_Material = new Material(renderingShader);
...中略...
// GPUバッファをマテリアルに設定
renderingShader_Material.SetBuffer("posBuffer", posBuffer);
次にSurfaceShader.shader側の頂点シェーダーをみると
	StructuredBuffer<double3> posBuffer;

	// 頂点シェーダ
	VSOut vert(uint id : SV_VertexID)
	{
		// idを元に、ドットの情報を取得
VSOut output; output.pos = mul(UNITY_MATRIX_VP, float4(float3(posBuffer[id]), 1)); return output; }

posBufferはdouble型なのでfloatにキャストしている。生の座標値だとまだ描画に使えないのでUNITY_MATRIX_VPでうまい具合に描画用の座標に変換している(正直よくわかってない・・・)。
ここでidには例によって0から始まる通し番号が入っている。これを使って各ドットの座標が格納されているposBufferにアクセスしている。

そして
C#側
    void OnRenderObject()
    {
        // レンダリングを開始
        renderingShader_Material.SetPass(0);
        // n個のオブジェクトをレンダリング
        Graphics.DrawProceduralNow(MeshTopology.Points, particleNum);
    }
ここのparticleNumで指定している数だけドットが表示されるというわけだ。

MeshTopology.Lines

lorenzline

MeshTopology.PointsをMeshTopology.Linesとかに変えてみると、こんな感じに違った雰囲気の絵が作れるのでいろいろ遊んでみると面白い。他にもQuadsなどもあり本格的なポリゴン描画にも使えそうである。

補足:GPUで擬似乱数の生成

実は今回、パーティクル座標の初期値をGPUで生成している。しかし大量の初期座標を重複なく作るのは意外と大変。Unity C#ではRandom.Random()などあるが、Compute Shaderではもちろん使えない。
そこでスレッドidから擬似ランダムな値を生成するコードを書いている。
//簡易ランダム
uint wang_hash(uint seed)
{
	seed = (seed ^ 61) ^ (seed >> 16);
	seed *= 9;
	seed = seed ^ (seed >> 4);
	seed *= 0x27d4eb2d;
	seed = seed ^ (seed >> 15);
	return seed;
}

//ランダム位置に初期値を置く
[numthreads(64, 1, 1)]
void BufInit(uint id : SV_DispatchThreadID)
{
	double3 pos;
	uint r = wang_hash(id * 1847483629);
	pos.x = 0.005 * (r % 1024) - 1.5;
	r /= 1024;
	pos.y = 0.005 * (r % 1024) - 3.5;
	r /= 1024;
	pos.z = 0.005 * (r % 1024) + 23.0;
	posBuffer[id] = pos;
}

この部分。まずスレッドidに1847483629という適当な数字をかけ、それをwang_hashという関数にかけている。たったこれだけ。
参考にしたのはQuick And Easy GPU Random Numbers In D3D11

この記事にもあるようwang_hashは乱数生成としては簡便な割には優秀らしい。実際TESTU01という乱数検定ソフトでも、「スレッドidに1847483629をかけてwang_hashに入れる」方式は最低限のテストをパスできた。



wikipediaから抜粋すると
『TESTU01には(10個のテストから成る)"Small Crush", (96個のテストからなる)"Crush"そして(106個のテストよりなる)"Big Crush"など数種のテスト・スイートが含まれる。』


もちろんXOR SHIFTにはかなわないが、XOR SHIFTは愚直に実装するとまったくの逐次計算になるのでGPUにはあまり向かない。一応、調べてみるとGPUに実装している人もちらほらいるのでできなくはないようだが、関数一つ書けば終わる程度のものとは比較にならないほど大変そうだ。


※本ブログで紹介しているソースコードはすべてCC0 パブリックドメインです。
ご自由にお使いください