«地理院地図 3D を Oculus Rift で見てみる 最新 髪の毛 (1)»

床井研究室

※このブログは遅くとも 2027 年 3 月に管理者の定年退職により閉鎖します (移転先は管理者本人共々模索中)

■ 2014年09月13日 [OpenGL] Oculus Rift でリアルタイムボリュームレンダリング

2021年10月30日 13:48更新

プログラミング

先日,久しぶりにあった知人に「現役プログラマなんですね」と言われてしまったけど,プログラミングは自分にとって表現の手段や考える道具みたいなものなので,こういう仕事をしている限りプログラミングと縁が切れることはないと思います.こういう歳になるとプログラミングは学生さんなんかに任せるべきっていう話になるのかもしれないけど,自分でプログラミングをしてないと発想そのものが出てきません.それにプログラミングをすればするほど自分の未熟さを思い知らされるので,やっぱり終わらない気がします.いつかはやめる時が来るのかもしれないけど,それまでは学生さんたちとガチでタメを張る感じでプログラミングしていたいと思います.

3D テクスチャ

ずいぶん前 (今見たら 8 年前!) に 3D テクスチャを使ったプログラムを書いた*1んですけど,その時に「これでボリュームレンダリングもできるじゃん」って思いました.でも,その方法は (当時使っていた) 教科書にも載ってる良く知られた方法だったので,結局実装しないまま時間が経ってしまいました.先日来 Oculus Rift で遊んでいて,Oculus で雲の中に突入したらどんな感じかなと思って作ってみました.

オリジナルのアイデアもわずかばかり入ってるし,これを出発点としてまた色々考えることもできそうなので,ちゃんとした業績につなげる前にソースを公開するってのは職業的にどうなのって思わないわけでもないんですけど,ほとんどは既に知られた内容ですし,大したプログラムでもないので,気にしないことにします.

ボリュームデータ

ボリュームデータは CPU 側で作りました (main.cpp の makeVolume() 関数).ボリュームのサイズが 32 × 32 × 32 だとデータの作成にあまり時間はかかりませんが,128 × 128 × 128 だと結構待たされます.GLSL にも noise 関数はあるので,それを使ってシェーダで作るべきなんでしょうけど,ここでは前述の以前作ったプログラムを流用しました.条件コンパイルってのは今風ではないんですけど,こらえてつかあさい.

    // 作業用メモリ
    std::vector<GLubyte> texture;
    
    // ノイズ関数を初期化する
    const Noise3 noise(5, 5, 5);
    
    // ノイズ関数を使ってテクスチャを作る
    for (GLint k = 0; k < depth; ++k)
    {
      const double z(double(k) / double(depth - 1));
      
      for (GLint j = 0; j < height; ++j)
      {
        const double y(double(j) / double(height - 1));
        
        for (GLint i = 0; i < width; ++i)
        {
          const double x(double(i) / double(width - 1));
          
#if VOLUMETYPE == CHECKER
          texture.push_back(((i >> 4) + (j >> 4) + (k >> 4)) & 1 ? 0 : 255);
#elif VOLUMETYPE == NOISE
          texture.push_back(GLubyte(noise.noise(x, y, z) * 255.0));
#elif VOLUMETYPE == PERLIN
          texture.push_back(GLubyte(noise.perlin(x, y, z, 4, 0.5) * 255.0));
#elif VOLUMETYPE == TURBULENCE
          texture.push_back(GLubyte(noise.turbulence(x, y, z, 4, 0.5) * 255.0));
#elif VOLUMETYPE == SPHERE
          const double px(2.0 * x - 1.0), py(2.0 * y - 1.0), pz(2.0 * z - 1.0);
          texture.push_back(GLubyte(255.0 - sqrt(px * px + py * py + pz * pz) * 127.5));
#else
          texture.push_back(255);
#endif
        }
      }
    }

VOLUMETYPE == TURBULENCE だと,こんな具合のボリュームテクスチャができます.

ボリュームデータ

これを 3D テクスチャに突っ込みます*2.ポイントは GL_CLAMP_TO_BORDER を指定して,境界色のアルファ値を 0 にしておくあたりです.

    // テクスチャオブジェクトを作成して結合する
    GLuint tex;
    glGenTextures(1, &tex);
    glBindTexture(GL_TEXTURE_3D, tex);
    
    // テクスチャを割り当てる
    glTexImage3D(GL_TEXTURE_3D, 0, GL_R8, width, height, depth, 0,
      GL_RED, GL_UNSIGNED_BYTE, &texture[0]);
    
    // テクスチャの拡大・縮小に線形補間を用いる
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    
    // テクスチャからはみ出た部分には境界色を用いる
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
    
    // テクスチャの境界色を設定する (ボリュームの外には何もない)
    static const GLfloat black[] = { 0.0f, 0.0f, 0.0f, 0.0f };
    glTexParameterfv(GL_TEXTURE_3D, GL_TEXTURE_BORDER_COLOR, black);

アルファブレンディング

このテクスチャをポリゴンにマッピングします.その際,ボリュームのセルの濃度をアルファ値にも設定しておきます.そしてアルファブレンディングを有効にすると,境界色のアルファ値を 0 にしているので,周囲の余計な部分が消えてボリュームテクスチャのところだけが見えるようになります.なお,アルファテストと違って,見えないところもレンダリングされています.

3D テクスクチャとしてマッピング

ポリゴンを複数枚並べるとボリューム全体を表示できそうです.アルファブレンディングを有効にすると,できそうな気がしてきます.

並べたポリゴンに 3D テクスチャをマッピング

スライスの描画

glDrawArraysInstanced() や glDrawElementsInstanced() を使えば,同じポリゴンを複数同時に描くことができます.ここでは,このポリゴン群をスライスと呼んでいます.

    // スライスの描画
    glBindVertexArray(slice);
    glDrawArraysInstanced(GL_TRIANGLE_FAN, 0, 4, slices);

同じところに複数のポリゴンを描いても仕方ないので,バーテックスシェーダで描画するインスタンスごとにポリゴンの位置をずらします (slice.vert).インスタンスは GLSL の組み込み変数 gl_InstanceID で識別できますから,これをもとにポリゴンの z 値を決定します.spacing はポリゴンの間隔で,1 / (ポリゴン数 - 1) です.gl_InstanceID に 0.5 を足しているのはテクスチャのサンプリング位置をセルの中心にするためで,ポリゴン群の中心 xy 平面上に移すために最後に 0.5 を引いています.

#version 150 core
#extension GL_ARB_explicit_attrib_location : enable
 
// テクスチャ座標の変換行列
uniform mat4 mt;
 
// モデルビュー変換行列
uniform mat4 mw;
 
// プロジェクション変換行列
uniform mat4 mp;
 
// スライスの間隔
uniform float spacing;
 
// [-0.5, 0.5] の正方形の頂点位置
layout (location = 0) in vec2 pv;
 
// スライスの頂点位置
out vec4 p;
 
// スライスのテクスチャ座標
out vec3 t;
 
void main()
{
  // スライスを gl_InstanceID でずらす
  p = vec4(pv, (float(gl_InstanceID) + 0.5) * spacing - 0.5, 1.0);
テクスチャ座標はスライスの頂点座標をもとに決定します.ボリュームの回転表示はスライスを回転させるのではなく,テクスチャ座標の方を回転します.回転によってボリュームがスライスからはみ出ないように,回転後のテクスチャ座標を √3 倍 (すなわちマッピングされるボリュームのサイズを 1 / √3 倍) します.
  // スライスのテクスチャ座標はスライスの中心を基準に √3 倍に拡大してから回転する
  t = (mat3(mt) * p.xyz) * 1.732 + 0.5;
  
  // 頂点位置を視点座標系に移す
  p = mw * p;
  
  // モデルビュープロジェクション変換をしてからラスタライザに送る
  gl_Position = mp * p;
}

しかし,アルファブレンディングを有効にしても,並べるポリゴンの枚数が多くなると,結局中身が見えなくなってしまいます.

ポリゴンが多いと中身が見えない

フラグメントの破棄

そこで,閾値を決めて,それを下回るアルファ値 (濃度) がマッピングされたフラグメントをフラグメントシェーダで破棄するようにします (slice.frag).

#version 150 core
#extension GL_ARB_explicit_attrib_location : enable
 
... (中略) ...
 
// テクスチャのサンプラ
uniform sampler3D volume;
 
// クリッピング座標系への変換行列
uniform mat4 mc;
 
// 閾値
uniform float threshold;
 
// スライスの表面上の位置
in vec4 p;
 
// テクスチャ座標
in vec3 t;
 
// フレームバッファに出力するデータ
layout (location = 0) out vec4 fc;
 
void main()
{
  // 元のボリュームの濃度と閾値の差
  float v = texture(volume, t).r - threshold;
  
  // 濃度が閾値以下ならフラグメントを捨てる
  if (v <= 0.0) discard;

すると,こんな風になります.

濃度が閾値以下のフラグメントを破棄

濃度勾配を求める

破棄しなかったフラグメントに関しては,そのフラグメントにマッピングされるセルの位置における濃度勾配を求めます.

  // 濃度の勾配を求める
  vec3 g = vec3(
    textureOffset(volume, t, ivec3(-1, 0, 0)).r - textureOffset(volume, t, ivec3(1, 0, 0)).r,
    textureOffset(volume, t, ivec3(0, -1, 0)).r - textureOffset(volume, t, ivec3(0, 1, 0)).r,
    textureOffset(volume, t, ivec3(0, 0, -1)).r - textureOffset(volume, t, ivec3(0, 0, 1)).r
  );

これを正規化したあと,0.5 倍して 0.5 を足して [0, 1] の範囲に直してフラグメントの色に使ってみます.

  // 勾配をそのままフラグメントの色に使う
  fc = vec4(normalize(g) * 0.5 + 0.5, v);

これは,こんな具合になります.この段階では意味はないけど,アルファブレンディングもしてみました.

濃度勾配を求める

陰影付け

濃度勾配はそのまま法線ベクトルとして使えるので,それを使って陰影付けを行います.

  vec3 l = normalize((pl * p.w - p * pl.w).xyz);  // 光線ベクトル
  vec3 n = normalize(g * mat3(mt));               // 法線ベクトル
  vec3 h = normalize(l - normalize(p.xyz));       // 中間ベクトル
  
  // 拡散反射光+環境光の反射光
  vec4 idiff = max(dot(n, l), 0.0) * kdiff * ldiff + kamb * lamb;
  
  // 鏡面反射光
  vec4 ispec = pow(max(dot(n, h), 0.0), kshi) * kspec * lspec;
  
  // フラグメントの色
  fc = vec4((idiff + ispec).rgb, v);
}

最終的には,こんな具合になります.キーボードの B のキーでアルファブレンディングの有効/無効を切り替えられます

陰影づけ

課題

すべてのスライスで勾配や陰影を求めているので,結構重いです.ボリュームデータが変化しないなら,勾配は事前に求めておくこともできます.これで結構速くなります*3.光源とボリュームデータ (およびスペキュラ*4を含めるなら視点) の位置関係が変化しないなら,さらに陰影も事前に計算しておくことができます.Forward Shading って言うんでしょうか.他にも色々考えてることはあるんですけど,まあ,私が思いつく程度のことだし,どうでもいいや.

勾配の事前計算

勾配を事前計算とフレームごとの処理時間の計測を追加しました (9 月 15 日).勾配を事前に計算するには slice.frag の記号定数 GRADIENT を 0 にしてください.

#define GRADIENT 0  // 勾配を事前計算しないなら 1

フレームごとの処理時間を表示するなら,config.h の記号定数 BENCHMARK を 1 にしてください.

// 経過時間を表示するなら 1
#define BENCHMARK     1

*1 この記事のコメントに OpenGL / GLSL には逆行列を求める方法が無いって書いてますけど,今は GLSL に inverse() っていう組み込み関数がありますね.

*2 ボリュームテクスチャを GLubyte (unsigned char) で作ってしまったので,テクスチャの internal format は GL_R8 (8 ビット 1 チャンネル) にしてます.

*3 最初は事前に勾配を求めていたのですが,一つのシェーダでできちゃうなぁと思ってまとめてしまいました.

*4 煙とか雲とかではあんまり考慮されることはないと思いますけど…

コメント(4) [コメントを投稿する]
還暦プログラマー 2015年01月24日 11:56

私もCGのプログラミングをしていて、ここの記事は読み物として重宝しています。<br>プログラミングは物書きと同じで年齢は関係ないと思っています。ミックジャガーより年上の酒井幸市先生がお手本です。<br>今後特にGPGPUの分野でご指導いただければ幸いです。

床井 2015年01月26日 15:27

還暦プログラマー様、コメントありがとうございます。<br>酒井幸市先生の本は私のところでも活用させていただいております。<br>こちらこそご指導いただければ幸いです。<br>よろしくお願いいたします。

しろ 2021年10月22日 11:03

お世話になっております.<br>ポリゴンを用いた可視化の部分について質問なのです.ポリゴンの面のなっている部分(フラグメントの濃度勾配を求めている画像でいう手前の面)の実行結果は納得がいきます.<br>しかし上下左右の面は複数のポリゴンの面ではなくエッジによって構成されているのでなぜ最終的に面になっているのか理解ができません.伝わりにくい質問になっていたら申し訳ありません.<br>お忙しいところ恐縮ですがよろしくお願いいたします.

とこ 2021年10月30日 13:48

しろさま、コメントありがとうございます。<br>確かに、上下左右の面は幾何学的にはエッジなのですが、ラスタライズしてしまえば画素の集まりにすぎません。<br>画素ごとに連続した陰影をつければ、その領域は面のように知覚されます。<br>でも、スライスの感覚が大きいと、向きによってはエッジの筋が見えます。<br>https://twitter.com/tokoik/status/1454306821909270528<br>スライスの方を視線に垂直のまま固定して、テクスチャ(座標)の方を回転させるべきだと思います。


編集 «地理院地図 3D を Oculus Rift で見てみる 最新 髪の毛 (1)»