GPUを使って、mrubyでのaobenchを更に速くした

おはようございます。先週は、次男のトビヒがひどくなり、奥さんが皮膚科に連れて行き、 その指示に従い、自分が、毎日次男の患部に薬をつけ、包帯を巻いていました。 それが、長男にもうつった疑惑があり、週の後半は長男もプール見学させてました。 @kjunichiです。

背景

aobenchの高速化にあたり、julialangでマルチプロセス化するも、パフォーマンス悪化、そこでmruby-aobench-nativeを作った。 折角なのでと、pthread使ってマルチスレッド化を試みるも、これまた、思うようなパフォーマンスが出なかった。 何とかしたい!

GPUのちからを借り、mruby-glslを作った。

mruby-glslとは

GLSLをmrubyのヒアドキュメントに書いて、これのレンダリング結果をバイナリ形式のppmデータとして 取得できるmrbgem。

そこそこの演算が必要になる画像データをGLSLによりGPUの高速な演算能力を利用して高速に取得可能となる。

ほんとに速いのか？

aobenchで比較

• ピュアrubyのaobench
• OSX

time bin/mruby benchmark/bm_ao_render.rb >ba.ppm

real    7m16.417s
user    7m10.905s
sys 0m2.536s

出力画像サイズを256x256に変更して実行。

• 先日のmruby-aobench-native
• OSX

time bin/mruby aon.rb > aon.ppm

real    0m2.991s
user    0m2.960s
sys 0m0.019s

• 今回のmruby-glsl
• OSX

aoglsl.rbの中身

vshader = <<'EOT'
#version 120

#ifdef GL_ES
precision mediump float;
#endif
vec4 qmul(vec4 q1, vec4 q2) {
        return vec4(
                dot(q1, q2.wzyx * vec4(1.0, 1.0, -1.0, 1.0)),
                dot(q1, q2.zwxy * vec4(-1.0, 1.0, 1.0, 1.0)),
                dot(q1, q2.yxwz * vec4(1.0, -1.0, 1.0, 1.0)),
                dot(q1, q2 * vec4(-1.0, -1.0, -1.0, 1.0)));
}
vec4 conj(vec4 q) {
        return q * vec4(-1.0, -1.0, -1.0, 1.0);
}
vec4 angle_axis(float rad, vec3 axis_f3) {
        float halfRad = rad * 0.5;
        vec3 f3 = normalize(axis_f3);
        vec4 result = f3.xyzz * sin(halfRad);
        result.w = cos(halfRad);
        return result;
}
vec3 rot3(vec4 q, vec3 v3) {
        vec4 f4 = v3.xyzz * vec4(1.0, 1.0, 1.0, 0.0);
        vec4 result = qmul(qmul(q, f4), conj(q));
        return result.xyz;
}
attribute vec3 aVertexPosition;
uniform vec2 angles;
varying vec3 org,dir;
void main() {
        vec3 center = vec3(0, 0, -3.5);
        vec4 q = qmul(angle_axis(angles.x, vec3(0.0, 1.0, 0.0)), angle_axis(angles.y, vec3(1.0, 0.0, 0.0)));
        gl_Position=vec4(aVertexPosition, 1);
        org=center + rot3(q, -center);
        dir=normalize(rot3(q, -vec3(-aVertexPosition.x*1.6,-aVertexPosition.y,1)));
}
EOT

fshader = <<'EOT'
#version 120
#ifdef GL_ES
precision mediump float;
#endif
// FragmentProgram
//
// porting GLSL by kioku based on syoyo's AS3 Ambient Occlusion
// [http://lucille.atso-net.jp/blog/?p=638]
varying vec3 org,dir;
struct Ray
{
vec3 org;
vec3 dir;
};
struct Sphere
{
vec3 center;
float radius;
};
struct Plane
{
vec3 p;
vec3 n;
};
struct Intersection
{
  float t;
  vec3 p;     // hit point
  vec3 n;     // normal
  int hit;
};
void shpere_intersect(Sphere s,  Ray ray, inout Intersection isect)
{
  // rs = ray.org - sphere.center
  vec3 rs = ray.org - s.center;
  float B = dot(rs, ray.dir);
  float C = dot(rs, rs) - (s.radius * s.radius);
  float D = B * B - C;
  if (D > 0.0)
  {
    float t = -B - sqrt(D);
    if ( (t > 0.0) && (t < isect.t) )
    {
        isect.t = t;
        isect.hit = 1;
        // calculate normal.
        vec3 p = vec3(ray.org.x + ray.dir.x * t,
            ray.org.y + ray.dir.y * t,
            ray.org.z + ray.dir.z * t);
        vec3 n = p - s.center;
        n = normalize(n);
        isect.n = n;
        isect.p = p;
      }
  }
}
void plane_intersect(Plane pl, Ray ray, inout Intersection isect)
{
// d = -(p . n)
// t = -(ray.org . n + d) / (ray.dir . n)
float d = -dot(pl.p, pl.n);
float v = dot(ray.dir, pl.n);
if (abs(v) < 1.0e-6)
  return; // the plane is parallel to the ray.
  float t = -(dot(ray.org, pl.n) + d) / v;
  if ( (t > 0.0) && (t < isect.t) )
  {
    isect.hit = 1;
    isect.t   = t;
    isect.n   = pl.n;
    vec3 p = vec3(ray.org.x + t * ray.dir.x,
          ray.org.y + t * ray.dir.y,
          ray.org.z + t * ray.dir.z);
    isect.p = p;
  }
}
Sphere sphere[3];
Plane plane;
void Intersect(Ray r, inout Intersection i)
{
for (int c = 0; c < 3; c++)
{
  shpere_intersect(sphere[c], r, i);
}
plane_intersect(plane, r, i);
}
void orthoBasis(out vec3 basis[3], vec3 n)
{
basis[2] = vec3(n.x, n.y, n.z);
basis[1] = vec3(0.0, 0.0, 0.0);
if ((n.x < 0.6) && (n.x > -0.6))
  basis[1].x = 1.0;
else if ((n.y < 0.6) && (n.y > -0.6))
  basis[1].y = 1.0;
else if ((n.z < 0.6) && (n.z > -0.6))
  basis[1].z = 1.0;
else
  basis[1].x = 1.0;
basis[0] = cross(basis[1], basis[2]);
basis[0] = normalize(basis[0]);
basis[1] = cross(basis[2], basis[0]);
basis[1] = normalize(basis[1]);
}
int seed = 0;
float random()
{
seed = int(mod(float(seed)*1364.0+626.0, 509.0));
return float(seed)/509.0;
}
vec3 computeAO(inout Intersection isect)
{
  const int ntheta = 8;
  const int nphi   = 8;
  const float eps  = 0.0001;
  // Slightly move ray org towards ray dir to avoid numerical probrem.
  vec3 p = vec3(isect.p.x + eps * isect.n.x,
                isect.p.y + eps * isect.n.y,
                isect.p.z + eps * isect.n.z);
  // Calculate orthogonal basis.
  vec3 basis[3];
  orthoBasis(basis, isect.n);
  float occlusion = 0.0;
  for (int j = 0; j < ntheta; j++)
  {
    for (int i = 0; i < nphi; i++)
    {
        // Pick a random ray direction with importance sampling.
        // p = cos(theta) / 3.141592
        float r = random();
        float phi = 2.0 * 3.141592 * random();
        vec3 ref;
        ref.x = cos(phi) * sqrt(1.0 - r);
        ref.y = sin(phi) * sqrt(1.0 - r);
        ref.z = sqrt(r);
        // local -> global
        vec3 rray;
        rray.x = ref.x * basis[0].x + ref.y * basis[1].x + ref.z * basis[2].x;
        rray.y = ref.x * basis[0].y + ref.y * basis[1].y + ref.z * basis[2].y;
        rray.z = ref.x * basis[0].z + ref.y * basis[1].z + ref.z * basis[2].z;
        vec3 raydir = vec3(rray.x, rray.y, rray.z);
        Ray ray;
        ray.org = p;
        ray.dir = raydir;
        Intersection occIsect;
        occIsect.hit = 0;
        occIsect.t = 1.0e+30;
        occIsect.n = occIsect.p = vec3(0, 0, 0);
        Intersect(ray, occIsect);
        if (occIsect.hit != 0)
          occlusion += 1.0;
      }
  }
// [0.0, 1.0]
occlusion = (float(ntheta * nphi) - occlusion) / float(ntheta * nphi);
return vec3(occlusion, occlusion, occlusion);
}
void main()
{
sphere[0].center = vec3(-2.0, 0.0, -3.5);
sphere[0].radius = 0.5;
sphere[1].center = vec3(-0.5, 0.0, -3.0);
sphere[1].radius = 0.5;
sphere[2].center = vec3(1.0, 0.0, -2.2);
sphere[2].radius = 0.5;
plane.p = vec3(0,-0.5, 0);
plane.n = vec3(0, 1.0, 0);

Intersection i;
i.hit = 0;
i.t = 1.0e+30;
i.n = i.p = vec3(0, 0, 0);

Ray r;
r.org = org;
r.dir = normalize(dir);
seed = int(mod(dir.x * dir.y * 4525434.0, 65536.0));

vec4 col = vec4(0,0,0,1);
Intersect(r, i);
if (i.hit != 0)
{
  col.rgb = computeAO(i);
}

gl_FragColor = col;
}
EOT


glsl = Glsl.new
glsl.attachVertexShader vshader
glsl.attachFragmentShader fshader

puts glsl.render

time bin/mruby aoglsl.rb >aoglsl.ppm

real    0m1.011s
user    0m0.074s
sys 0m0.056s

CPU実行に比べて、速くなった。

mruby-glslを通して学べたこと

インスタンス変数の取得

rubyで定義したクラスの中のインスタンス変数にCからアクセスする。

mrb_intern_litでインスタンス変数名をmrubyのシンボル型に変換して、 これをmrb_iv_getに指定することで、インスタンス変数が取得出来た。

mrb_value shader;
shader = mrb_iv_get(mrb, obj, mrb_intern_lit(mrb, "@vertexShader"));

Cの文字列とmrubyの文字列の相互変換

mrubyからC

int len;
const char *source;

len = RSTRING_LEN(shader);
source = RSTRING_PTR(shader);

Cからmruby

int size;
const char *buf="hogefuga";
mrb_value mrbStr;

mrbStr = mrb_str_new(mrb, buf, stelen(buf));

snprintfの返り値がヤバい

書き込んだバイト数ではなく、書き込んだとしたら 何バイト使うかを返す。

ググって出てきたＩＢＭサンプルが、paiza.ioで動かず、調査してわかった。 どうも、昔のsnprintfとC99あたりで規格ができてからのsnprintfは 挙動が違うものがあるのかも？

glReadPixelsでOpenGLの描画内容の取得

遅いと評判のglReadPixelsは大昔に使ったことがあるハズなので、 一旦これで実装してみた。

WindowsでMSVC++環境下での外部ライブラリを参照するmrbgemsの作りか

OSX以外でも試したく、Windowsにも対応してみた。

今回はWindows環境で、glfw3+Openglな描画環境でうごくコードが扱えた。

• mruby-glsl/mrbgem.rake at master · kjunichi/mruby-glsl · GitHub

に成果が詰まってる。

OpenGLの新しめのバージョンをglfwで利用するには

まだ、これはmrbgemには取り込んでいないが、

glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);

特に

glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);

が重要、これがないとOpenGLの初期化で失敗した。

参考資料

GLSLのコードは以下を流用させてもらいました。 - github.com

• aobench_mini.pdf - Google ドライブ

関連記事

• aobenchをjuliaでやってみた
• mrubyでaobenchを動かし、GIFアニメを出力できるようにした
• mruby-webcamをヘッドレス対応した

10年前の記事

• 昨日の検索ワード

2年前の記事

• JavaScriptをDockerを使って黒い画面に表示する