WebGLで綺麗な線香花火を作る

夏が完全に終わってしまいましたね。今年は花火をできなかった人も多かったことと思います。
そんなこともあってWebGL(three.js)を使って線香花火を作ってみました。
WebGL線香花火デモページ
ソースコード(github)

はじめに

写真や、スローモーションで撮影された動画などが参考になるのでまずは見てみましょう。
画像検索:線香花火の写真 参考にした動画

玉から飛び出した火花が分裂している、そして火花の動きが速すぎて人の目には点ではなく線に見えるのがわかると思います。
また、写真ではボケが遠近感を出していて美しいので再現したいですね。
この記事では、WebGLによる線香花火の描画について、個人的に線香花火の大事だと思う3つの要素

• 滑らかな軌跡
• 綺麗な枝分かれ
• 美しいボケ

それぞれについて解説します

滑らかな軌跡: 火花の軌跡の方程式

シンプルな自由落下ではなく、空気抵抗も入れた式を使うと、火花が急減速する様子を再現できたり風の効果を入れやすくて良いです。

運動方程式
F = ma
重力加速度
g = { x: 0, y: 0, z: -9.8 }
火花に加わる力は空気抵抗(簡単のために速度に比例することにする)と重力
F = -kv + mg
風も考慮に入れた場合
F = -k(v-wind) + mg
空気抵抗の係数kは粒子の断面積に比例することにします(質量の2/3乗に比例)
k = pow(m,2/3)
火花の粒子が小さくなるほど空気抵抗も小さくなりますがそれ以上に質量が小さくなるので減速しやすくなります。

この運動方程式から、火花の粒子の動きを

for(略){
  dt = 1 / 60 # 1フレームの時間
  position += velocity * dt
  velocity += (-k * (velocity - wind) / m + gravity) * dt
}

のように繰り返しで計算することもできますが、この運動方程式を解いて一発で時刻tの位置・速度を計算する方が便利です。

function velocityAt(v0, t) {
  return 時刻tでの速度を一発で出す計算式
}
function positionAt(p0, v0, t) {
  return 時刻tでの位置を一発で出す計算式
}

dv/dt = -k/m*(v-wind)+g
↓ 解くと
v(t) = (v0 - wind - g / k) * exp(-kt / m) + wind + m * g / k
↓ 積分すると
p(t) = p0 - m / k * (v0 - wind - m * g / k) * (exp(-kt / m) - 1) + (wind + m * g / k) * t

よく見る自由落下の式(放物線)と違う形をしていそうに見えますが、kを0に近づけると放物線とほぼ一致します。

綺麗な枝分かれ: 火花の分裂と時間経過

火花の速度がフレームレートに比べて非常に速いので、点ではなく1フレームの間に移動する軌跡を描くことにします。

// 1フレーム(時刻0〜時刻timeまで)の間に火花が通る軌跡を計算・描画
function renderAndProceedParticle(position, velocity, duration, time) {
  if (duration >= time) {
    // duration(寿命)>=timeの場合、時刻timeまで時間を進める
    render(position, velocity, dt)
    const nextPos = positionAt(position, velocity, time)
    const nextVel = velocityAt(position, velocity, time)
    particlesInNextFrame.push([nextPos, nextVel, duration - time])
  } else {
    // duration(寿命もしくは分裂)<timeの場合
    // 時刻をduration進める
    // 寿命なら終了、分裂なら分裂後の新しい火花を生成して時間を(time - duration)進める
    render(position, velocity, duration)
    const endPos = positionAt(position, velocity, duration)
    const endVel = velocityAt(position, velocity, duration)
    if (N個に分裂する場合) {
      N個の新しい火花(endPos, endVel).forEach(([pos, vel, dur]) => {
        renderAndProceedParticle(pos, vel, dur, time - duration)
      })
    }
  }
}

分裂のルール

分裂数はランダムに・分裂後の火花の質量もばらけるようにすると良いです。

また、作用反作用を無視した不自然な分裂は目立つので、少なくとも運動量保存の法則は満たしておくと良いです。

// m[i]: 分裂後の火花粒子の質量
// v[i]: 分裂後の火花粒子の速度(分裂前との相対速度)
v[i] = random() // このままだと sum(m[i]*v[i])=0 にならない
v[i] -= sum(m[i] * v[i]) / sum(m[i]) // これで総運動量が0になる

美しいボケ: 被写界深度

理想的な(中学の理科の教科書に載っているような)円形レンズのカメラでは、点光源のボケはベタ塗りの円になります。(画像検索:玉ボケ)
また、レンズの形を星型・五角形・六角形にするとボケもその形になります。
では火花の軌跡のボケをどうレンダリングするか考えてみましょう。
まず単純なケースとして、ボケの半径が一定・火花の明るさも一定の場合を考えます。

時刻t=0〜1の間を細かく区切り、各瞬間の火花のボケを全部描画して加算合成すればできそうですね。

ではある1つのピクセルに注目してみましょう。

このピクセルは、時刻t=0.23~0.52の間ボケの円内にあるので、明るさは0.52-0.23=0.29になります。
WebGLを使ってこのピクセルに0.23や0.52を書き込むには円柱(Sylinder)を使います。
円柱を軌跡に沿って、さらに断面が常にカメラを向くように変形させ、色として時刻を書き込んでみると


表面をレンダリングしたら0.52が、裏面をレンダリングしたら0.23が、注目しているピクセルに書き込まれます。
レンダリングしたい値は0.52と0.23の差なので、裏面のときは符号を逆転させ、浮動小数点テクスチャに対して加算描画すると、最終的に目的の値である0.29をピクセルに書き込むことができます。

これらの軌跡の計算・円柱の変形は全てVertexShaderでできます。

varying float vTime;
void main() {
  float t = position.z * time;
  float e = exp(-friction * t);
  vec3 wg = wind + gravity / friction;
  vec3 gpos = p + wg * t - (v - wg) * (e - 1.0) / friction;
  vec3 cpos = (viewMatrix * vec4(gpos, 1)).xyz;
  cpos += vec3(position.xy, 0) * mix(ra, rb, position.z) / cpos.z;
  gl_Position = projectionMatrix * vec4(cpos, 1);
  vTime = t;
}

uniform vec4 color;
varying float vTime;
void main() {
  gl_FragColor = (gl_FrontFacing ? vTime : -vTime) * color;
}

火花の明るさ・ボケ半径が一定でない場合は、色として書き込む値は時刻ではなく、

∫(時刻tでのボケの円の明るさ)dt = ∫light(t)*(1/distance(t)**2/radius(t)**2)dt
ただし
light(t): 時刻tでの火花の明るさ
distance(t): 時刻tでカメラからの距離
radius(t): 時刻tでのボケ半径

になります。(先ほどの明るさ・半径が一定の例ではこの積分が　定数*t になる)
このままだとうまく解けませんが、light(t)と(1/distance(t)**2/radius(t)**2)の部分をそれぞれ適当なtの多項式で近似すると積分結果もtの多項式になり、shaderで直接計算できるようになります。

// light(t)をtの2次関数で表現・(1/distance(t)**2/radius(t)**2)部分をtの1次関数で近似・distanceは面倒なので無視した場合
// (b0+b1*t+b2*tt) * ((1-t/time)/ra/ra+t/time/rb/rb)
float sb = t * (brightness0 + t * (0.5 * brightness1 + t * brightness2 / 3.0));
float sbt = t * t * (0.5 * brightness0 + t * (brightness1 / 3.0 + 0.25 * t * brightness2));
vBrightnessSum = ((sb - sbt) / ra / ra + sbt / rb / rb);

加算描画で裏面と表面の差を計算する性質上、depthTestは使えません。depthTestの結果、表面と裏面のうちの片方が描画されない場合、白飛びしたり暗くなったりします。
なので、depthTestをせずに、背景・後ろ側の火花・線香花火の玉や棒部分・手前側の火花、の順番でレンダリングしましょう。
レンズの形を五角形などにしたい場合は、円柱の代わりに五角柱を使えばOKです。ポリゴン数の削減にもなります。

その他

パラメータ調整

現実の値(重力加速度、火花の初速)を使うなど、物理的に正しく計算するのも大切ですが、見た目がそれっぽくなるよう感覚で適当にチューニングするのも大事です。(そもそも厳密な計算はできずに色々な近似が入ってるわけだし)
火花の速度は不自然じゃない程度に実物よりも遅くする方が飛び散っている感が出ます。

仕上げ

背景は真っ暗でも良いんですが、ぼんやりとボケている画像を使うと雰囲気が出ます。
線香花火の棒部分と玉・背景を加える、風速に強弱をつける、燃え始めの炎を追加、時間経過に伴って玉の形や大きさを変える・火花の数や分裂の仕方を変化させるなどの演出を入れて完成です。