LoginSignup
29
21

More than 5 years have passed since last update.

@aa_debdeb(Atsushi Asakura)

WebGL2でGray-Scott反応拡散系シミュレーション

Posted at

WebGL2でGray-Scott反応拡散系シミュレーションを実装しました。
reaction-diffusion.gif
実際に動いているものはこちらで確認することができます。
https://aadebdeb.github.io/WebGL_GrayScott_ReactionDiffusion_Simulation/
ソースコードはGithubに置いておきました。
aadebdeb/WebGL_GrayScott_ReactionDiffusion_Simulation: Gray-Scott Reaction Diffusion Simulation in WebGL

Gray-Scott反応拡散系はライフゲームと同じような格子ベースのシミュレーションです。各格子がその位置の値を保持するので、Gray-Scott反応拡散系では2成分$u$, $v$の値を保持することになります。WebGLを用いたシミュレーションではテクスチャをシミュレーション空間全体、テクスチャの各ピクセルを1つの格子とします。

今回はcanvas全体でシミュレーションを行うために、以下のようにcanvasと同じ大きさのテクスチャを作成しています。Gray-Scott反応拡散系では各格子が$u$, $v$の2つの実数値を保持するため、texImage2Dのformat引数をgl.RGに、type引数をgl.FLOATにしています。。テクスチャのx要素が$u$、y要素が$v$に対応しています。internalformat引数は32ビット精度のgl.RG32Fにしています。float値を書き込みできるようにEXT_color_buffer_float拡張を有効化しています。

function createFramebuffer(gl, sizeX, sizeY) {
  const framebuffer = gl.createFramebuffer();
  gl.bindFramebuffer(gl.FRAMEBUFFER, framebuffer);
  const texture = gl.createTexture();
  gl.bindTexture(gl.TEXTURE_2D, texture);
  gl.texImage2D(gl.TEXTURE_2D, 0, gl.RG32F, sizeX, sizeY, 0, gl.RG, gl.FLOAT, null);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
  gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, texture, 0);
  gl.bindFramebuffer(gl.FRAMEBUFFER, null);
  gl.bindTexture(gl.TEXTURE_2D, null);
  return {
    framebuffer: framebuffer,
    texture: texture
  };
}

const canvas = document.getElementById('canvas');
const gl = canvas.getContext('webgl2');
gl.getExtension('EXT_color_buffer_float');

let stateFbObjR = createFramebuffer(gl, canvas.width, canvas.height);
let stateFbObjW = createFramebuffer(gl, canvas.width, canvas.height);
const swapFramebuffer = function() {
  const tmp = stateFbObjR;
  stateFbObjR = stateFbObjW;
  stateFbObjW = tmp;
};

後述するように、シミュレーション自体はシェーダー用いてテクスチャにオフスクリーンレンダリングすることで行います。その際に1つ前の時間ステップの状態を保持しているテクスチャを参照する必要があります。オフスクリーンレンダリング先のテクスチャに現在参照しているテクスチャは使用できないので、フレームバッファを読み込み用(R=Read)、書き込み用(W=Write)の2つをそれぞれ作成して、swapFramebufferメソッドで入れ替えできるようにしています。

作成したフレームバッファにシミュレーションの初期状態を書き込みます。今回は、canvasの中心から一定距離以内はランダムな値、それ以外は$(u, v)=(0, 0)$となるようにしてます。書き込み後、swapFramebufferメソッドでフレームバッファを入れ替えておきます。

const initialize = function() {
  gl.bindFramebuffer(gl.FRAMEBUFFER, stateFbObjW.framebuffer);
  gl.useProgram(initializeProgram);
  gl.uniform2f(initializeUniforms['u_resolution'], canvas.width, canvas.height);
  gl.uniform2f(initializeUniforms['u_randomSeed'], Math.random() * 1000.0, Math.random() * 1000.0);
  renderToFillScreen();
  gl.bindFramebuffer(gl.FRAMEBUFFER, null);
  swapFramebuffer();
};

#version 300 es

precision highp float;

out vec2 o_state;

uniform vec2 u_resolution;
uniform vec2 u_randomSeed;

float random(vec2 x){
  return fract(sin(dot(x,vec2(12.9898, 78.233))) * 43758.5453);
}

void main(void) {
  vec2 st = (2.0 * gl_FragCoord.xy - u_resolution) / min(u_resolution.x, u_resolution.y);
  if (length(st) < 0.3) {
    o_state = vec2(
      random(gl_FragCoord.xy * 0.15 + u_randomSeed + vec2(231.32, 171.92)),
      random(gl_FragCoord.xy * 0.21 + u_randomSeed + vec2(131.17, 319.23))
    );
  } else {
    o_state = vec2(0.0);
  }
}

初期化が完了したので、レンダリングループ内でシミュレーションの進行と描画を行います。Gray-Scott反応拡散系シミュレーションは時間ステップを大きくしすぎると発散してしまうため、固定時間ステップでシミュレーションを行います。レンダリングループ内のupdateメソッドでシミュレーションをtimeStep秒だけ進め、renderメソッドでレンダリングをおこなっています。

let simulationSeconds = 0.0;

let previousRealSeconds = performance.now() * 0.001;
const loop = function() {

  const currentRealSeconds = performance.now() * 0.001;
  const nextSimulationSeconds = simulationSeconds + parameters['time scale'] * Math.min(0.2, currentRealSeconds - previousRealSeconds);
  previousRealSeconds = currentRealSeconds;

  const timeStep = parameters['time step'];
  while(nextSimulationSeconds - simulationSeconds > timeStep) {
    update(timeStep);
    simulationSeconds += timeStep;
  }
  render();

  animationId = requestAnimationFrame(loop);
}
loop();

requestAnimationFrameのフレームレートは環境により異なるので、経過時間をもとに複数回updateメソッドを呼び出すようにしています。実時間でシミュレーションを行うとなかなか発展していかないのでparameters[time scale]で時間の進む速度を大きくできるようにしています。parameters[time scale]が大きくなるほど一度のloopメソッドの呼び出しで実行するupdateメソッドの回数が多くなり、結果として負荷が大きくなるので注意してください。

updateメソッドでは以下のようにシミュレーションの次の状態を計算しています。u_stateTextureからその格子と近傍の現在の状態を取得して、陽解法で次の状態を求めています。今回は上端と下端、左端と右端が連続しているという境界条件になるようにしています。Gray-Scott反応拡散系のアルゴリズムはこの記事を参考にそのまま実装しています。

const update = function(deltaTime) {
  gl.bindFramebuffer(gl.FRAMEBUFFER, stateFbObjW.framebuffer)
  gl.useProgram(updateProgram);
  gl.activeTexture(gl.TEXTURE0);
  gl.bindTexture(gl.TEXTURE_2D, stateFbObjR.texture);
  gl.uniform1i(updateUniforms['u_stateTexture'], 0);
  gl.uniform2f(updateUniforms['u_diffusion'], parameters['diffusion U'], parameters['diffusion V']);
  gl.uniform1f(updateUniforms['u_feed'], parameters['feed']);
  gl.uniform1f(updateUniforms['u_kill'], parameters['kill']);
  gl.uniform1f(updateUniforms['u_timeStep'], deltaTime);
  gl.uniform1f(updateUniforms['u_spaceStep'], parameters['space step']);
  renderToFillScreen();
  gl.bindFramebuffer(gl.FRAMEBUFFER, null);
  swapFramebuffer();
};

#version 300 es

precision highp float;

out vec2 o_state;

uniform sampler2D u_stateTexture;
uniform float u_timeStep;
uniform float u_spaceStep;
uniform vec2 u_diffusion;
uniform float u_feed;
uniform float u_kill;

void main(void) {
  ivec2 coord = ivec2(gl_FragCoord.xy);
  ivec2 stateTextureSize = textureSize(u_stateTexture, 0);

  vec2 state = texelFetch(u_stateTexture, coord, 0).xy;

  vec2 left = texelFetch(u_stateTexture, ivec2(coord.x != 0 ? coord.x - 1 : stateTextureSize.x - 1, coord.y), 0).xy;
  vec2 right = texelFetch(u_stateTexture, ivec2(coord.x != stateTextureSize.x - 1 ? coord.x + 1 : 0, coord.y), 0).xy;
  vec2 down = texelFetch(u_stateTexture, ivec2(coord.x, coord.y != 0 ? coord.y - 1 : stateTextureSize.y - 1), 0).xy;
  vec2 up = texelFetch(u_stateTexture, ivec2(coord.x, coord.y != stateTextureSize.y - 1 ? coord.y + 1 : 0), 0).xy;

  vec2 laplacian = (left + right + up + down - 4.0 * state) / (u_spaceStep * u_spaceStep);

  o_state = state + u_timeStep * (u_diffusion * laplacian + vec2(
    state.x * state.x * state.y - (u_feed + u_kill) * state.x,
    -state.x * state.x * state.y + u_feed * (1.0 - state.y)
  ));
}

renderメソッドのJavaScriptコードとフラグメントシェーダー以下のようになっており、ここでスクリーンへの描画を行います。フラグメントシェーダー内ではu_targetとu_renderingでどの値を描画するか、どのように描画するかをそれぞれ決定しています。2Dの場合は値をそのまま描画し、3Dとして描画する場合は、この記事を参考に法線を計算してシェーディングしています。

const render = function() {
  gl.useProgram(renderProgram);
  gl.activeTexture(gl.TEXTURE0);
  gl.bindTexture(gl.TEXTURE_2D, stateFbObjR.texture);
  gl.uniform1i(renderUniforms['u_stateTexture'], 0);
  gl.uniform1i(renderUniforms['u_target'], parameters['target']);
  gl.uniform1i(renderUniforms['u_rendering'], parameters['rendering']);
  gl.uniform1f(renderUniforms['u_spaceStep'], parameters['space step']);
  renderToFillScreen();
}

#version 300 es

precision highp float;

out vec4 o_color;

uniform sampler2D u_stateTexture;
uniform int u_target;
uniform int u_rendering;
uniform float u_spaceStep;

float getValue(ivec2 coord) {
  vec2 state = texelFetch(u_stateTexture, ivec2(coord), 0).xy;

  if (u_target == 0) {
    return state.x;
  } else if (u_target == 1) {
    return state.y;
  } else {
    return abs(state.x - state.y);
  }
}

vec3 render2d(ivec2 coord) {
  return vec3(getValue(coord));
}

vec3 lambert(vec3 color, vec3 normal, vec3 lightDir) {
  return color * max(dot(normal, lightDir), 0.0);
}

vec3 render3d(ivec2 coord) {
  ivec2 stateTextureSize = textureSize(u_stateTexture, 0);
  float state = getValue(coord);
  float left = getValue(ivec2(coord.x != 0 ? coord.x - 1 : stateTextureSize.x - 1, coord.y));
  float right = getValue(ivec2(coord.x != stateTextureSize.x - 1 ? coord.x + 1 : 0, coord.y));
  float down = getValue(ivec2(coord.x, coord.y != 0 ? coord.y - 1 : stateTextureSize.y - 1));
  float up = getValue(ivec2(coord.x, coord.y != stateTextureSize.y - 1 ? coord.y + 1 : 0));

  vec3 dx = vec3(2.0 * u_spaceStep, 0.0, (right - left) / (2.0 * u_spaceStep));
  vec3 dy = vec3(0.0, 2.0 * u_spaceStep, (up - down) / (2.0 * u_spaceStep));

  vec3 normal = mix(normalize(cross(dx, dy)), vec3(0.0, 0.0, 1.0), 0.5);

  vec3 color = vec3(0.0);
  color += lambert(vec3(0.8), normal, vec3(1.0, 1.0, 1.0));
  color += lambert(vec3(0.3), normal, vec3(-1.0, -1.0, 0.3));
  return color;
}

void main(void) {
  vec2 state = texelFetch(u_stateTexture, ivec2(gl_FragCoord.xy), 0).xy;

  if (u_rendering == 0) {
    o_color = vec4(render2d(ivec2(gl_FragCoord.xy)), 1.0);
  } else {
    o_color = vec4(render3d(ivec2(gl_FragCoord.xy)), 1.0);
  }
}

WebGLでGray-Scott反応拡散系シミュレーションを行う方法について解説しました。GPUを使ってシミュレーションを進めることで、格子をピクセル単位で設定してもリアルタイムに計算することができました。

値をハイトマップと見なして法線を計算することで3D的にシェーディングする方法は、Polygon Lounge #1で見たもの(@kaiware007さんの作品?)がかっこよかったので真似しました。

NEORTに色と初期状態を変更したものを置いたのでそちらも見ていただけると幸いです。
https://neort.io/art/bilttjs3p9f9psc9oafg
reaction-diffusion2.gif

29
21
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
29
21