Ｂスプラインで折線近似
折線（黒線）をＢスプライン曲線（赤線）で近似する
サンプルプログラムを置きます。
http://g0307.hp.infoseek.co.jp/bspline.zip
Ｂスプラインの詳しい定義は以下サイトのように色々なところで
解説されています。
http://markun.cs.shinshu-u.ac.jp/learn/cg/cg5/index2.html
ただしこれは制御点を与えて描画するというもので
任意のデータをＢスプラインで近似する制御点を
求めるという類の解説はあまりないと思います。
今回、3次Ｂスプラインでノットが等間隔で近似対象折線が
周期的である場合に簡単にＢスプラインで近似する方法を述べます。
実際使う分には3次・等間隔で十分な場合が多いはずです。
3次周期Ｂスプラインは以下のように、一つの基底関数を
3つの区間に場合分けして定義されます。
任意の曲線を重み付き基底関数和で最小自乗で近似する理論は
以下の行列問題を解くことになります。これは別にＢスプライン
に限った話ではありません。
上式を解けば終了です。
Ｂスプラインの場合、定義台が狭いので基底関数同士の内積で
与えられる計量行列Ｂの要素はほとんどがゼロになります。
3次Ｂスプラインの場合は５重対角行列で、周期的な場合は
行列の右上、左下に折り返し成分がきます。
また、周期的な場合は基底関数はどれも同じ形のため
基底関数同士の内積値は以下の?〜?のパターンの組み合わせで
計算できます。
次にベクトルＦの計算方法です。
近似対象がＢスプラインのノットと同じ間隔で並んでいる
折線データの場合、内積値は以下３つのパターンの組み合わせで
計算できます。
これでＢ，Ｆともに求まったのであとはＢの逆行列を
求めればいいのですが、制御点列が多くなってくると
計算時間がかかりそうだし、行列演算処理を作るのが
めんどくさいですよね（実際には、周期５重対角行列を
少ないメモリで効率よく解くアルゴリズムがあるかもしれませんが・・・）
そこで、厳密な近似結果はいらないからなんとなく
近似できてそうな結果が欲しい、コードもあまり長く書きたくない
という場合の計算が以下です。
単純に反復計算で近似したＢスプラインと対象折線の二乗誤差の積分値を
減らすように更新をかけます。上のサンプルプログラムはこれで
解いてますが３２点の折線データで数回の反復で収束します。
プログラムは以下となります。
float y[128]; //近似したい折線データ float ba[128];//求めるＢスプライン制御点の値 float d[128]; //ワーク変数 // 反復更新(5回) for(it = 0;it < 5;it++){ for(i = 0;i < N;i++){ d[i] = 0; }; for(i = 0;i < N;i++){ i1 = i; i2 = (i+1)%N; i0 = (i - 1 + N) % N; // 差分更新 d[i2] += ( 0.0083333333 * ba[i0] + 0.1083333333 * ba[i1] + 0.05 * ba[i2] - 0.04166666667 * y[i0] - 0.125 * y[i1]); d[i1] += ( 0.1083333333 * ba[i0] + 0.45 * ba[i1] + 0.1083333333 * ba[i2] - 0.3333333333 * y[i0] - 0.3333333333 * y[i1] ); d[i2] += ( 0.05 * ba[i0] + 0.1083333333 * ba[i1] + 0.0083333333 * ba[i2] - 0.125 * y[i0] - 0.04166666667 * y[i1] ); }; for(i = 0;i < N;i++){ ba[i] -= d[i]; };
