閑古鳥

オールドプログラマの日記。プログラミングとか病気(透析)の話とか。

B-spline

ラグランジュより仕事で使っているものを理解しておいた方が良いかと思うので次は B スプライン補間をしてみます。

式を見て初めて B スプラインは外挿を行えない理由がわかりました。ただ、ノットベクトルの方がよくわかりません。任意の区間のデータを補間する場合、ノットベクトルには X データをごにょごにょして突っ込めば良いのでしょうか? とりあえずそれらしい結果にはなったのですが……。

NPlot でプロットするために、今日は C# で書きました。 C++ と大差ありませんが、 C# はよく知らないのでちょっといい加減。

>

class BSpline
{
public void Input(double[] x, double[] y, int N, int K)
{
x_.Clear();
y_.Clear();

x_.AddRange(x);
y_.AddRange(y);

degree_ = K;
int p = degree_ + x_.Count;
int m = degree_ - 1;
for (int i = 0; i <= p; ++i)
{
if (i < m)
{
knot_.Add(0);
}
else if (i >= x_.Count)
{
knot_.Add(x_[x_.Count - 1]);
}
else
{
knot_.Add(x_[i - m]);
}
}
}

public double Calculate(double x)
{
double result = 0.0;
for (int i = 0, size = y_.Count; i < size; ++i)
{
result += y_[i] * Basis(i, degree_, x);
}
return result;
}


private double Basis(int j, int n, double t)
{
if(n == 1) {
return knot_[j] <= t && t < knot_[j + 1] ? 1.0 : 0.0;
}
else {
double a = 0.0;
if ((knot_[j + n - 1] - knot_[j]) != 0.0)
{
a = (t - knot_[j]) / (knot_[j + n - 1] - knot_[j]) * Basis(j, n - 1, t);
}
double b = 0.0;
if ((knot_[j + n] - knot_[j + 1]) != 0.0)
{
b = (knot_[j + n] - t) / (knot_[j + n] - knot_[j + 1]) * Basis(j + 1, n - 1, t);
}
return a + b;
}
}

private int degree_;
private List<double> x_ = new List<double>();
private List<double> y_ = new List<double>();
private List<double> knot_ = new List<double>();
}

f:id:wata_d:20060502012314j:image

で、書いてみました。赤い線が B スプラインで、黒い線が昨日のラグランジュ補間です。ちょっと色が見えにくいですが……。

B スプラインは途中の実測点を通らないのでこんな線になります。