プログラミング
先日、久しぶりにあった知人に「現役プログラマなんですね」と言われてしまったけど、プログラミングは自分にとって表現の手段や考える道具みたいなものなので、こういう仕事をしている限りプログラミングと縁が切れることはないと思います。こういう歳になるとプログラミングは学生さんなんかに任せるべきっていう話になるのかもしれないけど、自分でプログラミングをしてないと発想そのものが出てきません。それにプログラミングをすればするほど自分の未熟さを思い知らされるので、やっぱり終わらない気がします。いつかはやめる時が来るのかもしれないけど、それまでは学生さんたちとガチでタメを張る感じでプログラミングしていたいと思います。
Oculus でリアルタイムボリュームレンダリング from Kohe Tokoi on Vimeo.
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 にしているので、周囲の余計な部分が消えてボリュームテクスチャのところだけが見えるようになります。なお、アルファテストと違って、見えないところもレンダリングされています。
ポリゴンを複数枚並べるとボリューム全体を表示できそうです。アルファブレンディングを有効にすると、できそうな気がしてきます。
スライスの描画
glDrawArraysInstanced() や glDrawElementsInstanced() を使えば、同じポリゴンを複数同時に描くことができます。ここでは、このポリゴン群をスライスと呼んでいます。
// スライスの描画
glBindVertexArray(slice);
glDrawArraysInstanced(GL_TRIANGLE_FAN, 0, 4, slices);
同じところに複数のポリゴンを描いても仕方ないので、バーテックスシェーダで描画するインスタンスごとにポリゴンの位置をずらします (slice.vert)。インスタンスは GLSL の組み込み変数 gl_InstanceID で識別できますから、これをもとにポリゴンの z 値を決定します。spacing はポリゴンの間隔で、$1 / (\text{ポリゴン数} - 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);
テクスチャ座標はスライスの頂点座標をもとに決定します。ボリュームの回転表示はスライスを回転させるのではなく、テクスチャ座標の方を回転します。回転によってボリュームがスライスからはみ出ないように、回転後のテクスチャ座標を $\sqrt{3}$ 倍 (すなわちマッピングされるボリュームのサイズを $1 / \sqrt{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 を足して $\left[0, 1\right]$ の範囲に直してフラグメントの色に使ってみます。
// 勾配をそのままフラグメントの色に使う
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