non vorrei lavorare

2020年度からの小学校プログラミング教育の必修化を親として迎えるブロガーの書く、子供との日常

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

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

背景

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

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

f:id:kjw_junichi:20160809072340p:plain

mruby-glslとは

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

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

ほんとに速いのか?

aobenchで比較

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の返り値がヤバい

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

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

glReadPixelsでOpenGLの描画内容の取得

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

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

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

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

に成果が詰まってる。

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

関連記事

10年前の記事

2年前の記事