GLSLのレイマーチングで地形を描画する

GLSLのレイマーチングで地形を描画する方法です。

地形を描画するには、レイをちょっとずつ進めてレイの位置が地形の高さよりも低くなったら、そこをレイと地形の交点とします。

地形を生成するための関数は次のようになり、平面座標を入力として高さを出力します。
ここでは、異なる周波数のサイン波を複数重ねたもので地形を作っています。

float terrain(in vec2 p) {
	p *= 0.2;
	float height = 0.0;
	float amp = 1.0;
	for (int i = 0; i < 10; i++) {
		height += amp * sin(p.x) * sin(p.y);
		amp *= 0.5;
		p *= 2.07;
	}
	return height * 5.0; // 適当な値で高さをスケール
}

レイを進めて行く処理は次のようになります。
ループ処理の中で、先ほどのterrain関数から高さを取得して、レイの位置と地形の距離を計算します。
もし、距離が一定以下になった場合は、交差したとみなしてループ処理を中断します。そうでない場合は、距離をもとにレイを進めます。
ここで注意しないといけないのは、距離はレイの進行方向と地形の距離ではなく、あくまでy軸方法との距離だということです。そのため、断崖絶壁のような高さが急激に変わる地形ではノイズが目立つようになると思います。
ここでは、１未満の値を距離にかけて慎重にレイを進めるようにしています。

const float tmax = 100.0;
float raymarch(in vec3 origin, in vec3 direction) {
	float t = 0.0;
	for (int i = 0; i < 64; i++) {
		vec3  p = origin + t * direction;
		float height = terrain(p.xz);
		float distance = p.y - height;
		if(distance < 0.02 || t > tmax) {break;}
		t += 0.2 * distance; // 適当な値で距離を小さくする
	}
	return t;
}

法線は以下のようにレイマーチングでよく使われる数値差分で求めることができます。

vec3 normal(in vec2 p) {
	float epsilon = 0.02;
	return normalize(vec3(
		terrain(p + vec2(epsilon, 0.0)) - terrain(p - vec2(epsilon,0.0)),
		2.0 * epsilon,
		terrain(p +vec2(0.0, epsilon)) - terrain(p - vec2(0.0, epsilon))
	));
}

法線をもとにカラーリングを行うとこうなります。

全コードと動いているところは以下から確認できます。
http://glslsandbox.com/e#44275.0

ここからは、細かいテクニック的な話をします。

raymarch関数の中で、レイと地形との交差判定に以下の式を用いました。

if(distance < 0.02 || t > tmax) {break;}

これを以下のように、distanceとの比較にtを用いることで、近いところほど正確に、遠いところは粗い交差判定を行うことができます。

if(distance < 0.001 * t || t > tmax) {break;}

同様に、法線の計算でも以下のように交差点のtを渡すことで、近いところほど細かい計算を行うようにできます。

vec3 normal(in vec2 p, float t) {
	float epsilon = 0.001 * t;
	return normalize(vec3(
		terrain(p + vec2(epsilon, 0.0)) - terrain(p - vec2(epsilon,0.0)),
		2.0 * epsilon,
		terrain(p +vec2(0.0, epsilon)) - terrain(p - vec2(0.0, epsilon))
	));
}

他のテクニックとして、次のように交点を一つ前のステップの値と補間することでより正確なtの値を求める方法があります。

const float tmax = 100.0;
float raymarch(in vec3 origin, in vec3 direction) {
	float t = 0.0;
	float step = 0.0;
	float lastdistance = 0.0;
	for (int i = 0; i < 64; i++) {
		vec3  p = origin + t * direction;
		float height = terrain(p.xz);
		float distance = p.y - height;
		if(distance < 0.02) {
			t  += step * distance / (lastdistance - distance);
			break;
		}
		if (t > tmax) {break;}
		step = 0.2 * distance;
		t += step;
		lastdistance = distance;
	}
	return t;
}

詳しくは、次のサイトで解説されています。

地形のライティングは以下のリンク先の内容が参考になります。
http://www.iquilezles.org/www/articles/outdoorslighting/outdoorslighting.htm