標高データを2次元スプライン補間

前回で取り上げた地形のハイトマップ（標高データ）は、解像度が約100mで
あり、人間スケールのオブジェクトを置くには少々粗い。そこで、頂点を補完することを考えてみた。

単なる線形補間だと、曲率の高い地点で不自然な「折れ曲がり」が見えてしまうので、スムーズな補間を行いたい。この要求を満足するものに自然スプライン補間があった。

理論的なところは全くフォローできていないが、下のページなどを参考にして２次元スプライン補間を実装することができた。

http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal1/node16.html


これは補間を行いたい区間の両端と、その一つ外側のデータ点のみがあれば補間できるため、動的LODによってポリゴンを削減した部分まで計算に入れる必要が無く、比較的高速にオンメモリで計算できるという利点がある。ただし、地形のような莫大な数の点を補間するとなると計算のオーバーヘッドも無視できなくなってくるので、オフライン計算やキャッシュが有効であろう。

1枚目の画像はSRTM3をデータ源としてレンダリングした富士山頂付近をキャプチャしたものの補間前と後の比較である。画像では補間は16分割して行ったが、理論上はいくつにでも分割できる(厳密に言うと、ピラミッドバッファのサイズから上限が決まってくる)。

なお、実際の富士山頂にも登ってみたが、2枚目の画像のような光景だった。

リアルの地形データからなめらかな地形を描画したいときには、有効な方法ではなかろうか。

続きにスプライン補間のソースを掲示する。
 

/* 3次のスプライン補間、y座標を求める関数
   a: 4点のy座標
   x: 補間点のx座標
   返り値: 補間点のy座標 */

static double spline_cubic(const double a[4], double x){
	int i, j;
	double alpha[4], l[4], mu[4], z[4];
	double b[4], c[4], d[4];
	for(i = 1; i < 3; i++)
		alpha[i] = 3. * (a[i+1] - a[i]) - 3. * (a[i] - a[i-1]);
	l[0] = 1.;
	mu[0] = 0.;
	z[0] = 0.;
	for(i = 1; i < 3; i++){
		l[i] = 4. - mu[i-1];
		mu[i] = 1. / l[i];
		z[i] = (alpha[i] - z[i-1]) / l[i];
	}
	l[3] = 1.;
	z[3] = 0.;
	c[3] = 0.;
	for(j = 2; 0 <= j; j--){
		c[j] = z[j] - mu[j] * c[j+1];
		b[j] = a[j+1] - a[j] - (c[j+1] + 2. * c[j]) / 3.;
		d[j] = (c[j+1] - c[j]) / 3.;
	}
	return a[1] + b[1] * x + c[1] * x * x + d[1] * x * x * x;
}

/* 3次のスプライン補間、補間曲線の方程式を求める関数
   a: 4点のy座標
   f: 3次関数の係数
      y=a+bx+cx^2+dx^3 のとき
      f[0] == a
      f[1] == b
      f[2] == c
      f[3] == d
 */

static void spline_cubic_f(const double a[4], double f[4]){
	int i, j;
	double alpha[4], l[4], mu[4], z[4];
	double b[4], c[4], d[4];
	for(i = 1; i < 3; i++)
		alpha[i] = 3. * (a[i+1] - a[i]) - 3. * (a[i] - a[i-1]);
	l[0] = 1.;
	mu[0] = 0.;
	z[0] = 0.;
	for(i = 1; i < 3; i++){
		l[i] = 4. - mu[i-1];
		mu[i] = 1. / l[i];
		z[i] = (alpha[i] - z[i-1]) / l[i];
	}
	l[3] = 1.;
	z[3] = 0.;
	c[3] = 0.;
	for(j = 2; 0 <= j; j--){
		c[j] = z[j] - mu[j] * c[j+1];
		b[j] = a[j+1] - a[j] - (c[j+1] + 2. * c[j]) / 3.;
		d[j] = (c[j+1] - c[j]) / 3.;
	}
	f[0] = a[1];
	f[1] = b[1];
	f[2] = c[1];
	f[3] = d[1];
}