ElixirからOpenGLを使って3D空間に描画をする

この記事はElixir Advent Calendar 2020の13日目の投稿です。
昨日12日目の投稿は、「書評：プログラミングElixir第2版」でした！

はじめに

こんにちは。普段はライブ配信サービスのバックエンドエンジニアをしている @Sadalsuud です。

今回は、ElixirでCreative Codingみたいなことができるのではないかと思い立ってやってみたという内容で、こういう感じの描画結果をリアルタイムで得ようというものです。

ElixirについてもwxwidgetもOpenGLについても、甘い部分についてはお目こぼしいただければと思っています。内容はあくまで私の解釈です。リポジトリはここで公開しています。

前提/この記事で扱う範囲

Elixirからwxwidgetを通じてOpenGLを叩くのですが、基本的な描画を行うところまではすでにチュートリアルを作成してくださっている方がいらっしゃいます。ですので、本記事ではこのチュートリアルが意図していたことの解説から始まり、さらに冒頭でお見せした描画結果を表示するまでの追加の手順をざっくり見ていきたいと思います。

Getting started with OpenGL in Elixir

なお、上記チュートリアルを進めるにあたってハマった2点

• :wx.demoでGUIを起動しようとするとCould not find wxe_driver.so と言われる
• WSLから描画結果を表示させるための VcXsrv Windows X Server を導入する。

についてはこちらに対処法がありました。

また筆者は業務でElixirを使用したことはなくほぼ雰囲気で書いています。もし間違いやより適した書き方など、心当たりのある方がいらっしゃいましたらご一報いただければ大変うれしく思います。

描画にはそれなりに負荷がかかるようなので私の環境を記載しておきます。
- OS Windows 10 Pro WLS2 Ubuntu-18.04
- CPU Intel(R) Core(TM) i9-10900K
- GPU GeForce RTX 3080

Elixirからポリゴンを描画するまで

Getting started with OpenGL in Elixir
このチュートリアルを終えると、表示されたウィンドウに三角形のポリゴン板が表示されています。

たったこれだけなのですが、1loopごとに毎回ポリゴンの頂点が定義され描画を繰り繰り返しています。
全体の流れはこうなると思われます。

1. 最初に一度だけstart_link/0が呼ばれ、GUI/OpenGLの初期化が行われたのち、ループ間で受け渡される{frame, %{canvas: canvas, timer: timer}}というタプルが初期化されます。その中身が描画ごとに受け渡されていきます。

2. 描画の前にはhandle_info(:update, state)が呼ばれます。start_link/0の返り値として生成されたタプルもstateとして引数に渡されています。この中でさらにrender/1にstateが渡されます。

3. renderはdrawを呼びます。draw/0は描画にまつわる処理が書か書かれており、最初に現在描画に使われていないほうのバッファをクリアして、座標系をリセットします。そのあと座標系を並行移動させた後、ポリゴンの頂点を3か所定義しています。beginに渡しているのはsrc/gl_const.erlから呼び出しているGL_TRIANGLESというOpenGLの頂点をどの順番で結んでいくか定義する部分です。

4. drawを抜けたあとは現在表示中のバッファと、先ほどまで操作したバッファをswapして新しい描画結果を表示させます。

5. 次の描画が始まります。描画の前にはhandle_info(:update, state)が呼ばれます........

このような処理の流れで1 loopごとに三角形が描画されています。

:wxや :gl はerlangのモジュールをElixirから呼ぶときのもので、ドキュメントを見ることでそれぞれ使い方のヒントがあります
:gl -> https://erlang.org/doc/man/gl.html
:wx -> https://erlang.org/doc/man/wx.html

ちなみに、実装中何度もstart_linkを手打ちするのは大変なので.iex.exsにElixirOpengl.start_link()と入れておくと、$ iex -S mix したときに自動で実行開始してくれるようになって楽でした。

1loopごとに描画されていることを確かめる

現状ではloopが回っているとはいえ、描画ごとに変化がないのでウィンドウの中は静止画のように止まっています。
そこで、initの最後の行を次のように変更してみます。

lib/elixir_opengl.ex

  def init(config) do
    wx = :wx.new(config)
   # 中略...

    {frame, %{canvas: canvas, timer: timer, count: 0}} # 変更！
  end

今、stateに新たにcountというkeyを持たせ0をセットします。
そして次にhandle_info(:update, state)を次のように変更します。

lib/elixir_opengl.ex

  def handle_info(:update, state) do
    :wx.batch(fn -> render(state) end)
    next_state = Map.merge(state,%{count: state.count + 1})
    {:noreply, next_state}
  end

drawの中を変更し

lib/elixir_opengl.ex

  defp draw(state) do #stateを渡すようにする

    :gl.clear(Bitwise.bor(:gl_const.gl_color_buffer_bit, :gl_const.gl_depth_buffer_bit))
    :gl.loadIdentity()
    :gl.translatef(0.0, 0.0, -6.0) # 変更 x方向への移動を0にする
    :gl.rotatef(state.count, 0.0, 0.0, 1.0) # 追加 ベクトル x: 0.0 y: 0.0 z: 1.0 を軸にcountの量だけ回転
    :gl.'begin'(:gl_const.gl_triangles)
    :gl.vertex3f(0.0, 1.0, 0.0)
    :gl.vertex3f(-1.0, -1.0, 0.0)
    :gl.vertex3f(1.0, -1.0, 0.0)
    :gl.'end'()
    :ok
  end

こうすることで、1loopごとにstateの中にあるcountが1づつ増加していきます。
この増加が回転量となって三角形が回ります。

描画を行う関数を作ってみる

次は円を描いてみます。
円を描く関数を定義して、その関数をdrawの中で呼ぶことで描画を行います。

次のcircleという関数は次のようなタプルを引数にとり円を描画します。
{半径,頂点数,{赤,緑,青}}

lib/elixir_opengl.ex

  defp circle(radius, resolution, {r, g, b}) do # 新しく定義
    deg_to_rad = :math.pi() / 180
    :gl.begin(:gl_const.gl_polygon()) # 後でsrc/gl_const.erl に追加する
    :gl.color3f(r, g, b)

    Enum.map(
      0..(resolution + 1),
      fn x ->
        :gl.vertex3f(
          :math.cos(x * (360 / resolution) * deg_to_rad) * radius,
          :math.sin(x * (360 / resolution) * deg_to_rad) * radius,
          0.0
        )
      end
    )
    :gl.end()
  end

この関数でやっていることは、三角関数を使って半径radiusの円の軌跡上の点をresolution+1数だけ繰り返し取得しては頂点として登録している感じです。頂点数を上げていけばきれいな円が描画されます。+1しないと始点と終点がうまく閉じなかったのでそうしています。
ついてに色もつけています。

途中で:gl_const.gl_polygon()を使っているのですが、この部分は頂点を繋ぐ順番をOpenGLに伝えるもので、現状まだこんな値は定義されていないので、src/gl_const.erlに追記をします。

src/gl_const.erl

# 任意の場所に追加
gl_polygon() ->
  ?GL_POLYGON.

そしてdrawを書き換えcircleを呼ぶと晴れて円が描画されます。

lib/elixir_opengl.ex

  defp draw(state) do
    :gl.clear(Bitwise.bor(:gl_const.gl_color_buffer_bit, :gl_const.gl_depth_buffer_bit))
    :gl.loadIdentity()
    :gl.translatef(0.0, 0.0, -6.0) 
    :gl.rotatef(state.count, 0.0, 0.0, 1.0)
    circle(1,12,{0.0,0.2,0.7}) #{半径,頂点数,{赤,緑,青}}
    :ok
  end

座標と速度を保存する

あるフレームで描画された図形を少し動かしたいとき、前のフレームの座標 + 動かしたい量　のように処理して少しずつ場所を移動させたいことがあります。

これを実現するには、次の処理に渡すstateに、現在の座標を保存してから渡してあげれば、次のフレームで前のフレームの位置を参照できそうです。

今回の例では、粒子は座標と速度を持ちます。それぞれ3次元のベクトルで表すことができるので、一つの粒子は position{x,y,z},velocity{vx,vy,vz}を持つことになります。

これをまとめた構造体 Particle を作ります。

lib/particle.ex

defmodule Particle do
  defstruct position: {0, 0, 0}, velocity: {0, 0, 0}

  def init_particles(count) do
    Enum.map(
      0..count,
      fn _ ->
        %Particle{
          position:
            {:random.uniform(100) / 100, :random.uniform(100) / 100, :random.uniform(100) / 100},
          velocity: {0, 0, 0}
        }
      end
    )
  end
end

init_particles/1はランダムな座標でcount個分のPerticle構造体のリストを作成してくれる関数です。
次に、実際に使って初期化をします。

lib/elixir_opengl.ex

  def init(config) do
    wx = :wx.new(config)
   # 中略...

    {
     frame, %{canvas: canvas,
     timer: timer,
     count: 0,
     particles: Particle.init_particles(100)}
    } 
  end

こうすることで100個のParticle構造体のリストを作成してstateに保存することができました。

とりあえず表示してみるときはこんな感じです。

lib/elixir_opengl.ex

  defp draw(state) do
    :gl.clear(Bitwise.bor(:gl_const.gl_color_buffer_bit, :gl_const.gl_depth_buffer_bit))
    :gl.loadIdentity()
    :gl.translatef(0.0, 0.0, -3.0)
    Enum.map(state.particles, fn particle -> particle.position end)
    |>Enum.map(fn pos -> circle(0.03, 32, pos, {0.1, 1, 1}) end)
    :ok
  end

実行すると

この状態では動きに変化はありません。次のフレームに対して現在と同じParticle構造体を受け渡しているので、常に同じ位置に表示されています。粒子が右上の固まっているのはランダムに生成した座標がすべて正の数だからこうなっています。

次は1ループごとに少しずつ位置をずらしてみます。

ベクトル場を作って座標と速度に変化を与える

次はベクトル場を作って、1つひとつの粒子に速度変化を与えていきます。
もちろんベクトル場じゃなくてもいいのですが、単調に加算するだけだと、動きがあまりきれいじゃないのでベクトル場を用いています。

個人的にCurlNoiseが大好きなので、それっぽいものを作って行こうと思い、 「Three.jsでCurl Noise」 をたくさん参考にさせていただき、好みのベクトル場を目指しました。

ただ、このノイズを作るためにSimplexNoizeの生成をElixir上で再現しなければならないのですが、この処理をElixirに書き換えるのがつらくて心が折れてしまい、差し当たって公開されていたパッケージ noise を使ってしまいました。しかし、このnoise使えはするもののgithubのページも削除されており、どういう仕組みで生成しているのか詳しく追えていません。

が、いろいろひっかき回してなんとなくかわいい動きになったので、全部ヨシとさせてください

mix.exs

  defp deps do
    [
      {:noise, "~> 0.0.2"},
    ]
  end

lib/vector_field.ex

defmodule VectorField do
  def affect_particle(particle) do
    field_vel = curl_noise(particle.position)

    new_velocity = add3d(particle.velocity, field_vel) |> multiply3d({0.92, 0.92, 0.92})
    new_position = add3d(particle.position, multiply3d(new_velocity, {0.01, 0.01, 0.01}))

    %Particle{position: new_position, velocity: new_velocity}
  end

  def curl_noise({x, y, z} = pos) do
    multiply3d(pos, {3, 3, 3}) |> simplex_noise_delta3d |> multiply3d({0.1, 0.1, 0.1})
  end

  defp simplex_noise3d({x, y, z} = pos) do
    x_noise = :noise_simplex.raw(2, multiply3d(pos, {0.5, 0.5, 0.5} |> add3d({30, 30, 30})))
    y_noise = :noise_simplex.raw(3, multiply3d(pos, {0.5, 0.5, 0.5} |> add3d({30, 30, 30})))
    z_noise = :noise_simplex.raw(4, multiply3d(pos, {0.5, 0.5, 0.5} |> add3d({30, 30, 30})))
    {x_noise, y_noise, z_noise}
  end

  defp simplex_noise_delta3d({x, y, z} = pos) do
    dlt = 0.0001
    a = simplex_noise3d(pos)
    b = simplex_noise3d({x + dlt, y + dlt, z + dlt})

    {(elem(a, 0) - elem(b, 0)) / dlt, (elem(a, 1) - elem(b, 1)) / dlt,
     (elem(a, 2) - elem(b, 2)) / dlt}
  end

  defp random_add(n) do
    :random.uniform(10) / 100 + n
  end

  def add3d({x1, y1, z1}, {x2, y2, z2}) do
    {x1 + x2, y1 + y2, z1 + z2}
  end

  def sub3d({x1, y1, z1}, {x2, y2, z2}) do
    {x1 - x2, y1 - y2, z1 - z2}
  end

  def multiply3d({x1, y1, z1}, {x2, y2, z2}) do
    {x1 * x2, y1 * y2, z1 * z2}
  end

  defp normalize3d({x, y, z}) do
    len = length3d({x, y, z})

    if len == 0 do
      {0, 0, 0}
    else
      {x / len, y / len, z / len}
    end
  end

  defp length3d({x, y, z}) do
    :math.sqrt(:math.pow(x, 2) + :math.pow(y, 2) + :math.pow(z, 2))
  end
end

affect_particle/1はParticle構造体を渡して、新しいParticle構造体を返す関数です。
handle_info(:update, state)を書き換え、state.paritciles に入っているPerticle構造体の一つひとつをaffect_particleに渡しベクトル場の影響を与えていきます。

lib/elixir_opengl.ex

  def handle_info(:update, state) do
    :wx.batch(fn -> render(state) end)
     new_particles = state.particles
       |>Enum.map(fn particle -> VectorField.affect_particle(particle) end)
    next_state = Map.merge(state,%{count: state.count + 1,particles: new_particles})
    {:noreply, next_state}
  end

depth.getしたあとに実行すると、フレームごとに粒子の位置が変化します。

それぞれの座標と座標の距離を計算する

ここから、近い点と点を線で繋ぐような表現を目指してみたいと思います。愚直に計算すると1つの点からすべての点への距離を計算しないといけないため、n × n の計算量になってしまうと思います。

ですが今回は点Aと点Bの距離と、逆の点Bと点Aの距離は同じものと扱えるので省略可能なはずです。また、点Aと点Aの距離は常にゼロだし、描画する必要もないので計算する必要もありません。そう考えると計算すればいい対応は、表の〇の部分だけになると考えられます。

この計算を再帰っぽく関数にして落とし込んでみるとこのような形になりました。
(正しい再帰ではない気もしたので「っぽい」と言いました。もっと上手な方法を知りたいです。)

lib/connect_line.ex

defmodule ConnectLines do
  @cutoff 10

  def get_edges([result | []]), do: Enum.reject(result, fn x -> is_nil(x) end)

  def get_edges([result | array]) do
    [a | next_array] = array

    near_points = Enum.map(next_array, fn b -> if distance3d(a, b) < @cutoff, do: {a, b}, else: nil end)
    next_result = result ++ near_points
    get_edges([next_result | next_array])
  end

  def distance3d({x1, y1, z1}, {x2, y2, z2}) do
    length3d({x1 - x2, y1 - y2, z1 - z2})
  end

  def length3d({x, y, z}) do
    :math.sqrt(:math.pow(x, 2) + :math.pow(y, 2) + :math.pow(z, 2)) |> abs
  end
end

このget_edges関数は、下記のように、particle構造体から取り出した座標の配列を渡すことで、再帰的に計算を行うことを想定しています。

    ConnectLines.get_edges( [[] | particle_positions])

このようにget_edage呼ぶと、一番最初は get_edges([result | array])にマッチし、先頭の要素a と、aを取り除いたあとの配列arrayに分かれて処理が進行します。
aに対してarrayの各要素bとそれぞれ対応とをり、distance3d/2を使って3次元の2つの座標間の距離を求めます。その値が @coutoff の値より小さければ {a,b} 大きければ nil が 返され、これらのリストが result に追加されたうえで、次のループが回ります。

これが再帰的に実行され最終的には get_edges([result | []]) にマッチし、今までnilだった要素を取り除いたうえで結果を返します。

こうすると下記のようなタプルのリストが取得できます。

[
  {{x1,y1,z1},{x2,y2,z2}},
  {{x3,y3,z3},{x4,y4,z4}},
  {{x5,y5,z5},{x6,y6,z6}},
  ...
]

これは、距離が近い点のリストで、このリストの各タプルの1つ目の座標と2つ目の座標をそれぞれ繋いでいけば、全体として近い点がつながったような表現ができるはずです。

そこでdrawの中で線を描画する処理を加えます。

lib/elixir_opengl.ex

defp draw(state) do
    :gl.clear(Bitwise.bor(:gl_const.gl_color_buffer_bit, :gl_const.gl_depth_buffer_bit))
    :gl.loadIdentity()
    :gl.translatef(0.0, 0.0, -100.0)

    particle_positions = Enum.map(state.particles, fn particle -> particle.position end)
    Enum.map(particle_positions,fn pos -> circle(0.5, 32, pos, {0.1, 1, 1}) end)

    ConnectLines.get_edges( [[] | particle_positions])
    |> Enum.map(fn x -> Tuple.to_list(x) end)
    |> List.flatten
    |> lines(1)

    :ok
  end

  defp lines(positions, size) do
    :gl.pointSize(size)
    :gl.begin(:gl_const.gl_lines())

    Enum.map(
      positions,
      fn {x, y, z} -> :gl.vertex3f(x, y, z) end
    )
    :gl.end()
  end

gl_const.erl

gl_lines() ->
  ?GL_LINES.

drawの中でつなぐべき2つの座標のリストをConnectLines.get_edges( [[] | particle_positions]) から取得して、それらをフラットな座標リストにしたうえでline/1を使って線を描画します。

するとこのような感じの結果が表示されます。

並列化してみる

【注意】ここでは恐らく雑な並列化をしています。また測定方法も再現性のある測定方法をしていません。「どうやら処理が分散していそうだ」というところまでを荒く見ているのでそのつもりでお願いします。

Particleの数が100くらいでは描画は遅延なく行われますが、この数をどんどん上げていくと計算量は爆発していくのでどんどん実行が遅くります。私の環境では500くらいではもう映像はカクカクでした。

簡易的にフレームレートを計る

そこで並列に計算をできるようにして、フレームレートがどれだけ変化するか試してみます。
今回は簡易的に、更新回数 / 実行開始から経過した秒数 とすることでだいたいの実行速度を見れるようにします。実行開始時のstateに uptime: :os.system_time(:millisecond) のようにunixtimeを保存しておいて、updateのたびにプログラム起動から何秒経過したか調べ、今までのloopした回数を割ってログに出すだけです。

exec_time = (:os.system_time(:millisecond) - state.uptime) / 1000
IO.inspect(state.count / exec_time)

処理の並列化

今まで一つひとつの粒子にベクトル場を適応していた場所はこのように並列化できそうです。

lib/elixir_opengl.ex

    # 並列
    new_particles =
      state.particles
      |> Flow.from_enumerable(stages: 18)
      |> Flow.map(fn particle -> VectorField.affect_particle(particle) end)
      |> Enum.to_list()

    # 非並列
    # new_particles = state.particles
    # |>Enum.map(fn particle -> VectorField.affect_particle(particle) end)

そして、計算量の多かった点と点の距離を計算する部分を並列にしてみます。

connect_lines.ex

 def get_edges_parallel(array) do
    indexed_array = Enum.with_index(array)

    near_points =
      indexed_array
      |> Flow.from_enumerable(stages: 18)
      |> Flow.map(fn n -> sort_positions(n, indexed_array) end)
      |> Enum.to_list()
      |> List.flatten()
      |> Enum.reject(fn x -> is_nil(x) end)
  end

  def sort_positions(target, array) do
    Enum.map(array, fn x ->
      if elem(x, 1) < elem(target, 1) && distance3d(elem(target, 0), elem(x, 0)) < @cutoff,
        do: {elem(target, 0), elem(x, 0)},
        else: nil
    end)
  end

lib/elixir_opengl.ex

    near_points = ConnectLines.get_edges_parallel(particle_positions) # 並列
    |> Enum.map(fn x -> Tuple.to_list(x) end)
    |> List.flatten
    |> lines(1)

    # ConnectLines.get_edges( [[] | particle_positions]) # 非並列
    # |> Enum.map(fn x -> Tuple.to_list(x) end)
    # |> List.flatten
    # |> lines(1)

このようにすることで実行の一部を並列に実行できそうです。

計測

粒子の数を500にして実行開始約30秒の時点での 更新回数/秒 を求めました。
実際に実行したのはリポジトリにあるコードなので、この記事のガイドで完成させたコードとは異なることを付記します。
また交絡なく比較するための特別な調整を行っていないので、あくまで断片的な情報です。

非並列時 -> 4.20009081277433
並列時 -> 13.601519571263822

となり、雑な計測なので再現性があるかはわかりませんが、並列実行時がたくさん画面更新をできているという結果になりました。
またその時にリソースモニタを見てみると、なんだか負荷が分散されているようにも見えます。（再現性は確認できていないので、雑感としてご覧ください。）

まとめ

これで少し回転や色の調整、ウィンドウの大きさなどをいじって完成です。

せっかく関数型言語なのだから時刻を引数にとって座標を出力する関数が書けれれば、わざわざstateで構造体を管理するようなこともしなくていいのかもしれません。並列化についても、考えなくFlowを使ってみた感じなのでまだまだ改善の余地がありそうです。
:wx_misc.getMousePosition() を使えばカーソルの位置情報も取れるのでインタラクティブな表現も可能だと思います。

また解決できないこととしてウィンドウを動かしたときにキャンバスの挙動がおかしいという問題があります。wxの問題かOpenGlの問題か、はたまた VcXsrv Windows X Serveの問題か切り分けられていません。もし心当たりある方は教えていただければとてもうれしく思います。

いろいろできてないことは多いですが、それでも絵を出すことと同じくらいElixirのコードでの実装を考えるのが楽しく感じましたし、パフォーマンスを出すための工夫ができることも興味の源泉になっています。

次は並行性をもっと生かせないかチャレンジしてみたいと思います。

明日は @mnishiguchi さんの「GenServerのプロセスをどう管理するか」です！お楽しみに！