Interpolation — 補間

概要

sangi の補間モジュールは、 既知の標本点 $(x_i, y_i)$ から間の値を再構成する関数群を提供する。 用途に応じて次のカテゴリに分かれる。

• 区分・多項式補間 — ラグランジュ / Newton 分割差分 / 線形・双線形 / 階段 (zero-order hold) / 重心ラグランジュ
• 3 次スプライン — 自然境界・クランプ境界の $C^2$ 連続な区分 3 次曲線
• 形状保存・局所補間 — PCHIP (単調保存)・修正 Akima・Catmull-Rom・滑らか補間・Floater-Hormann 有理補間・5 次 Hermite
• B-スプライン — Cox-de Boor 基底・クランプ節点・de Boor 評価・B-スプライン補間 / 回帰
• 三角補間 — 等間隔周期データの三角多項式補間
• スプラインフィッティング — 平滑化パラメータ付き最小二乗フィット (1D / パラメトリック / 2D テンソル積)
• 画像補間 — 最近傍・双 3 次 (Catmull-Rom カーネル)・Lanczos リサンプリング

テンプレートパラメータ $T$ は concepts::Field または concepts::OrderedField を要求する。 通常は double、 高精度が必要なら多倍長型 Float を渡せる。 すべて namespace sangi。

多くの関数は std::span<const T> を取る版と std::vector<T> を取る版の両オーバーロードを持つ。 本ページでは主に std::vector 版で説明する (span 版も引数の意味は同じ)。

ビルド

// 区分多項式・スプライン・形状保存・フィッティング・画像補間
#include <math/interpolation/Interpolation.hpp>

// De Boor 基底・cubicSpline・B-スプライン回帰
#include <math/interpolation/BSpline.hpp>

ヘッダオンリー。 リンクするライブラリは不要。 内部で線形系を解く関数 (B-スプライン補間・スプラインフィッティング・サーフェススプライン) は <math/linalg/solvers.hpp> に依存するが、 Interpolation.hpp が自動でインクルードする。

区分・多項式補間

関数詳細

lagrange_interpolation

template<concepts::Field T>
T lagrange_interpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T xi);

動作: $n$ 個の標本点を通る次数 $n-1$ のラグランジュ補間多項式 $L(x) = \sum_i y_i \prod_{j \ne i} \frac{x - x_j}{x_i - x_j}$ を点 $x_i$ で直接評価する。 各評価が $O(n^2)$。

パラメータ:

注意: サイズ不一致・空入力は std::invalid_argument を送出。 同じ標本点で多数の点を評価するなら、 重みを再利用できる barycentricLagrangeWeights / barycentricLagrangeEvaluate の方が速い。 高次の等間隔点では Runge 現象に注意 (floaterHormannWeights を参照)。

// 例: 3 点 (0,0),(1,1),(2,4) を通る放物線を x=1.5 で評価
std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
double v = lagrange_interpolation<double>(xs, ys, 1.5);
// 実行結果: v = 2.25  (= 1.5^2)

newton_interpolation_coefficients

template<concepts::Field T>
std::vector<T> newton_interpolation_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

動作: Newton の前進分割差分表を構築し、 補間多項式 $p(x) = c_0 + c_1 (x - x_0) + c_2 (x - x_0)(x - x_1) + \cdots$ の係数 $\{c_i\}$ を返す。 係数は newton_interpolation_evaluate に渡して評価する。

パラメータ:

用途: 標本点を固定して何度も評価する場合。 係数を一度求めれば各評価は Horner 法で $O(n)$。

newton_interpolation_evaluate

template<concepts::Field T>
T newton_interpolation_evaluate(
    const std::vector<T>& x,
    const std::vector<T>& coeffs,
    T xi);

動作: newton_interpolation_coefficients の返した分割差分係数を Horner 法で点 $x_i$ で評価する。

パラメータ:

std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
auto c = newton_interpolation_coefficients<double>(xs, ys);
double v = newton_interpolation_evaluate<double>(xs, c, 1.5);
// 実行結果: v = 2.25  (ラグランジュ補間と同じ多項式)

linear_interpolate

template<concepts::Field T>
T linear_interpolate(T x0, T y0, T x1, T y1, T xi);

動作: 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る直線上で $x_i$ の値を返す ($y_0 + t(y_1 - y_0)$, $t = (x_i - x_0)/(x_1 - x_0)$)。 $x_0 \approx x_1$ のときは縮退を避けるため 2 点の平均 $(y_0 + y_1)/2$ を返す。

パラメータ:

bilinear_interpolate

template<concepts::Field T>
T bilinear_interpolate(
    T x, T y,
    T x1, T y1, T x2, T y2,
    T q11, T q12, T q21, T q22);

動作: 矩形 $[x_1, x_2] \times [y_1, y_2]$ の 4 隅の値から点 $(x, y)$ の値を双線形補間する。 下辺・上辺で線形補間してから縦方向に再度線形補間する。

パラメータ:

注意: 評価点が矩形 $[x_1, x_2] \times [y_1, y_2]$ の外にある場合は std::invalid_argument を送出する。

stepInterpolation

template<concepts::OrderedField T>
T stepInterpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T xi);

動作: 階段補間 (zero-order hold)。 $x_i$ を含む区間 $[x_i, x_{i+1})$ の左端の値 $y_i$ をそのまま返す。 区間検索は二分探索で $O(\log n)$。

パラメータ:

注意: 範囲外は最も近い端点の値を返す ($x_i \le x_0$ で $y_0$、 $x_i \ge x_{n-1}$ で $y_{n-1}$)。

barycentricLagrangeWeights

template<concepts::Field T>
std::vector<T> barycentricLagrangeWeights(
    const std::vector<T>& x);

動作: 重心ラグランジュ補間の重み $w_j = 1 / \prod_{i \ne j} (x_j - x_i)$ を計算する。 重みは標本点だけで定まり、 評価点 $y$ には依存しないので一度計算すれば使い回せる。

パラメータ:

用途: 通常のラグランジュ補間は毎回 $O(n^2)$ だが、 重心形式は重み計算 $O(n^2)$ + 評価ごと $O(n)$。 同じ標本点で何度も評価する場合に有利。

barycentricLagrangeEvaluate

template<concepts::Field T>
T barycentricLagrangeEvaluate(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& weights,
    T xi);

動作: 第 2 種重心公式 $p(x_i) = \dfrac{\sum_j w_j y_j / (x_i - x_j)}{\sum_j w_j / (x_i - x_j)}$ で補間値を評価する。 $x_i$ がいずれかの標本点に一致する場合はその $y_j$ をそのまま返す。

パラメータ:

std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
auto w = barycentricLagrangeWeights<double>(xs);
double v = barycentricLagrangeEvaluate<double>(xs, ys, w, 1.5);
// 実行結果: v = 2.25

3 次スプライン

3 次スプラインは隣接区間を $C^2$ 連続で滑らかに接続する区分 3 次曲線である。 sangi は 2 つの入口を提供する。

• 係数配列を返す cubic_spline_coefficients + cubic_spline_evaluate (Interpolation.hpp)
• 係数と eval をまとめた結果構造体を返す cubicSpline (BSpline.hpp)

SplineBoundaryCondition

enum class SplineBoundaryCondition {
    Natural,   // 自然スプライン (両端で 2 階微分が 0)
    Clamped    // クランプスプライン (両端で 1 階微分値を指定)
};

cubic_spline_coefficients

template<concepts::Field T>
std::vector<std::array<T, 4>> cubic_spline_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    int boundary_condition = 0,
    const std::pair<T, T>& bc_values = { 0, 0 });

動作: 標本点を通る区分 3 次スプラインの係数を計算する。 各区間の係数は $\{a, b, c, d\}$ で、 区間 $i$ では $S_i(x) = a + b(x - x_i) + c(x - x_i)^2 + d(x - x_i)^3$。 内部で三重対角系を Thomas 法 (前進消去 + 後退代入) で解く。 区間数 $n-1$ 個の std::array<T, 4> を返す。

パラメータ:

注意: span 版は第 3 引数を SplineBoundaryCondition の列挙値で受ける。 サイズ不一致・点数 2 未満・$x$ が単調増加でない場合は std::invalid_argument を送出。 返した係数は cubic_spline_evaluate で評価する。

cubic_spline_evaluate

template<concepts::Field T>
T cubic_spline_evaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 4>>& coeffs,
    T xi);

動作: スプライン係数を点 $x_i$ で評価する。 $x_i$ を含む区間を線形探索し、 $S_i(x_i) = a + b\,\delta + c\,\delta^2 + d\,\delta^3$ ($\delta = x_i - x_i^{\text{left}}$) を返す。 範囲外は最終区間の定数項を返す。

パラメータ:

用途: cubic_spline_coefficients だけでなく、 pchipCoefficients・catmullRomCoefficients・modifiedAkimaCoefficients・smoothInterpolationCoefficients はいずれも同じ $\{a, b, c, d\}$ 形式の係数を返すため、 この関数で評価できる。

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 1.0, 0.0, 1.0};
auto coeffs = cubic_spline_coefficients<double>(xs, ys);  // 既定: 自然境界
double v = cubic_spline_evaluate<double>(xs, coeffs, 1.5);
// 実行結果: v ≒ 0.5 付近 (標本点 (1,1)-(2,0) の間を滑らかに通過)

cubicSpline

template<typename T>
struct CubicSplineResult {
    std::vector<T> a, b, c, d;  // 区間ごとの係数: S_i(x) = a + b(x-x_i) + c(x-x_i)^2 + d(x-x_i)^3
    std::vector<T> x;           // 節点
    T eval(T t) const;          // 補間値を評価
};

template<typename T>
CubicSplineResult<T> cubicSpline(
    const std::vector<T>& xs,
    const std::vector<T>& ys);

動作: BSpline.hpp が提供する自然 3 次スプライン補間 (両端で $S'' = 0$)。 係数と評価メソッド eval をまとめた CubicSplineResult<T> を返すので、 係数配列と $x$ を別々に持ち回る必要がない。

パラメータ:

返り値: CubicSplineResult<T>。 result.eval(t) で任意点の補間値を得る。

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 1.0, 4.0, 9.0};
auto spline = cubicSpline<double>(xs, ys);
double v = spline.eval(1.5);
// 実行結果: v ≒ 2.3 付近 (節点 (1,1)-(2,4) の間を滑らかに通過)

形状保存・局所補間

次の関数群は区間ごとの 3 次 (一部 5 次) 係数を返し、 cubic_spline_evaluate (5 次は quinticHermiteEvaluate) で評価する。 全体最適化を行わず局所的に傾きを決めるため、 標準 3 次スプラインより振動 (オーバーシュート) を抑えやすい。

関数詳細

pchipCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> pchipCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

動作: PCHIP (Piecewise Cubic Hermite Interpolating Polynomial)。 Fritsch-Carlson アルゴリズムで各節点の傾きを定め、 単調区間では単調性を保存してオーバーシュートを抑える。 内部点は調和平均、 端点は Bessel の片側 3 点公式で傾きを推定し、 単調性条件 $\alpha^2 + \beta^2 \le 9$ を満たすよう傾きを制限する。 区間ごとの $\{a, b, c, d\}$ を返し、 cubic_spline_evaluate で評価する。

パラメータ:

用途: 単調なデータ (累積分布・濃度プロファイル等) を滑らかに、 しかし非物理的なアンダー / オーバーシュートなしで補間したいとき。 点数 2 未満やサイズ不一致では空ベクトルを返す。

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 0.0, 1.0, 1.0};  // 段差状の単調データ
auto coeffs = pchipCoefficients<double>(xs, ys);
double v = cubic_spline_evaluate<double>(xs, coeffs, 1.5);
// 実行結果: 0 <= v <= 1 (単調保存ゆえ範囲外に飛び出さない)

modifiedAkimaCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> modifiedAkimaCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

動作: Akima の局所傾き推定に Boost.Math 準拠の修正を加えた版。 傾き差 $|\delta_{i+1} - \delta_i|$ の絶対値に微小量 $\varepsilon$ を足すことでゼロ除算を回避し、 平坦域・等間隔データでの不安定さを改善する。 区間ごとの $\{a, b, c, d\}$ を返し、 cubic_spline_evaluate で評価する。

パラメータ:

用途: 外れ値の影響を局所に閉じ込めたいデータ。 一点の異常値が遠くの区間まで波及しにくい (Akima の特徴)。 点数 2 未満やサイズ不一致では空ベクトルを返す。

catmullRomCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> catmullRomCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

動作: Catmull-Rom スプライン。 各節点の接線を両隣の節点から差分で求める $C^1$ 連続な補間スプライン。 端点は片側差分。 区間ごとの $\{a, b, c, d\}$ を返し、 cubic_spline_evaluate で評価する。

パラメータ:

用途: コンピュータグラフィックスでのキーフレーム補間・曲線描画。 点数 2 未満やサイズ不一致では空ベクトルを返す。

smoothInterpolationCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> smoothInterpolationCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

動作: 隣接区間の傾きを混合して接続を滑らかにする区分補間。 混合式 $A_{\text{mix}} = (|A_0| A_1 + |A_1| A_0)/(|A_0| + |A_1|)$ により、 3 点が一直線なら直線を再現し、 片側が平坦なら混合も 0、 山型 (傾きが逆符号) なら混合 0 となる。 先頭区間と末尾区間は端点傾きを拘束した 2 次、 中間区間は両端傾きを拘束した 3 次。 区間ごとの $\{a, b, c, d\}$ を返し、 cubic_spline_evaluate で評価する。

パラメータ:

注意: $x$ が狭義単調増加でない場合は std::invalid_argument を送出。

floaterHormannWeights

template<concepts::OrderedField T>
std::vector<T> floaterHormannWeights(
    const std::vector<T>& x,
    std::size_t d = 3);

動作: Floater-Hormann 重心有理補間の重みを計算する。 次数 $d$ の局所多項式をブレンドした有理補間で、 等間隔点でも Runge 現象を起こさず安定。 $d = 0$ で区分定数、 $d = n-1$ で多項式補間に対応する。

パラメータ:

用途: 高次の等間隔データを多項式補間したいが Runge 現象を避けたい場合。 重みは floaterHormannEvaluate に渡す。

floaterHormannEvaluate

template<concepts::OrderedField T>
T floaterHormannEvaluate(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& weights,
    T xi);

動作: 重心公式 $r(x_i) = \dfrac{\sum_k w_k y_k / (x_i - x_k)}{\sum_k w_k / (x_i - x_k)}$ で Floater-Hormann 有理補間を評価する。 評価点が節点に一致する (または極めて近い) 場合はその $y_k$ を返す。

パラメータ:

quinticHermiteCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 6>> quinticHermiteCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& dy,
    const std::vector<T>& d2y);

動作: 各節点で値 $f$・1 階微分 $f'$・2 階微分 $f''$ を指定する 5 次 Hermite 補間。 $C^2$ 連続。 区間ごとに 6 係数 $\{a, b, c, d, e, f\}$ を返し、 $S(t) = a + b\,t + c\,t^2 + d\,t^3 + e\,t^4 + f\,t^5$ ($t = x - x_i$) で評価する。

パラメータ:

用途: 曲率まで指定したい軌道・運動プロファイルの補間。 サイズ不一致・点数 2 未満では空ベクトルを返す。 評価には quinticHermiteEvaluate を使う (係数長が 6 なので cubic_spline_evaluate は不可)。

quinticHermiteEvaluate

template<concepts::OrderedField T>
T quinticHermiteEvaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 6>>& coeffs,
    T xi);

動作: quinticHermiteCoefficients の係数を点 $x_i$ で評価する。 区間は二分探索で特定し、 範囲外は端区間にクランプする。 Horner 法で 5 次多項式を評価。

パラメータ:

B-スプライン

B-スプラインは局所台 (compact support) を持つ基底関数の線形結合で曲線を表す方式。 sangi は 2 系統の基底評価を提供する。

• Interpolation.hpp の bsplineBasis(degree, i, x, knots) (Cox-de Boor 再帰、 OrderedField)・bsplineCoefficients・bsplineEvaluate
• BSpline.hpp の bsplineBasis(i, p, t, knots) (De Boor 再帰、 int 次数)・bsplineBasisAll・uniformKnots・deBoor・bsplineRegression

両ヘッダに同名の bsplineBasis があるが引数並びが異なる (前者は std::span<const T> の節点と std::size_t 次数、 後者は std::vector<T> の節点と int 次数)。

関数詳細

bsplineBasis (Interpolation.hpp)

template<concepts::OrderedField T>
T bsplineBasis(
    std::size_t degree,
    std::size_t i,
    T x,
    std::span<const T> knots);

動作: Cox-de Boor 再帰で B-スプライン基底関数 $B_{i,\text{degree}}(x)$ を評価する。 次数 0 では区間 $[\text{knot}_i, \text{knot}_{i+1})$ の指示関数 (最終区間のみ右端を含む)。

パラメータ:

bsplineCoefficients

template<concepts::OrderedField T>
std::pair<std::vector<T>, std::vector<T>> bsplineCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    std::size_t degree = 3);

動作: 標本点を補間する B-スプラインの節点ベクトルと制御点係数の組 {knots, coefficients} を返す。 両端に多重度 $\text{degree}+1$ のクランプ節点を置き、 内部節点は平均化節点で選ぶ。 コロケーション行列 $B_{j,\text{degree}}(x_i)$ を組んで LU 分解で線形系を解く。

パラメータ:

注意: サイズ不一致、 または点数が $\text{degree}+1$ 未満では std::invalid_argument を送出。 返した knots と coefficients は bsplineEvaluate に渡す。

bsplineEvaluate

template<concepts::OrderedField T>
T bsplineEvaluate(
    const std::vector<T>& knots,
    const std::vector<T>& coefficients,
    std::size_t degree,
    T x);

動作: 制御点係数と基底の線形結合 $\sum_i c_i B_{i,\text{degree}}(x)$ で B-スプライン補間値を評価する。

パラメータ:

std::vector<double> xs{0.0, 1.0, 2.0, 3.0, 4.0}, ys{0.0, 1.0, 0.0, 1.0, 0.0};
auto [knots, coeffs] = bsplineCoefficients<double>(xs, ys, 3);
double v = bsplineEvaluate<double>(knots, coeffs, 3, 2.5);
// 標本点を通る 3 次 B-スプライン上の値 (xs[2]=2, xs[3]=3 の間)

bsplineBasis (BSpline.hpp)

template<typename T>
T bsplineBasis(size_t i, int p, T t, const std::vector<T>& knots);

動作: De Boor 再帰で基底 $N_{i,p}(t)$ を評価する (Interpolation.hpp 版と同じ Cox-de Boor 漸化式だが、 次数を int で、 節点を std::vector で取る別オーバーロード)。

パラメータ:

bsplineBasisAll

template<typename T>
std::vector<T> bsplineBasisAll(int p, T t, const std::vector<T>& knots);

動作: 点 $t$ における全基底 $N_{0,p}(t), \ldots, N_{n-1,p}(t)$ ($n = \text{knots.size}() - p - 1$) を一括評価する。 $t$ が節点ベクトルの右端のときは最後の基底だけを 1 にする (右端の取りこぼし防止)。

パラメータ:

uniformKnots

template<typename T>
std::vector<T> uniformKnots(int n, int p, T a = T{0}, T b = T{1});

動作: 区間 $[a, b]$ にクランプ一様節点ベクトル (両端の多重度 $p+1$、 内部は等間隔) を生成する。 節点総数は $n + p + 1$。

パラメータ:

deBoor

template<typename T>
std::vector<T> deBoor(
    const std::vector<std::vector<T>>& controlPoints,
    const std::vector<T>& knots,
    int p,
    T t);

動作: De Boor アルゴリズムで B-スプライン曲線 $C(t) = \sum_i N_{i,p}(t) P_i$ 上の 1 点を評価する。 節点スパンを特定し、 三角表を再帰的に縮約して曲線上の点 (次元 $d$ のベクトル) を返す。

パラメータ:

返り値: 曲線上の点 (制御点と同じ次元 $d$ の std::vector<T>)。

bsplineRegression

template<typename T>
struct BSplineRegressionResult {
    std::vector<T> coefficients;  // B-スプライン係数
    std::vector<T> knots;         // 節点ベクトル
    int degree;                   // 次数
    T eval(T t) const;            // 予測値を評価
};

template<typename T>
BSplineRegressionResult<T> bsplineRegression(
    const std::vector<T>& x,
    const std::vector<T>& y,
    int nBasis = 10,
    int degree = 3);

動作: $n_{\text{basis}}$ 個の B-スプライン基底に対する最小二乗回帰。 一様節点で設計行列 $B$ を作り、 正規方程式 $B^\top B\, c = B^\top y$ を (微小正則化を加えて) Gauss 消去で解く。 補間 (全点を通る) ではなく平滑なフィット。 結果型の eval(t) で予測値を得る。

パラメータ:

用途: ノイズを含むデータの平滑な傾向曲線を引きたいとき。 基底本数 nBasis で滑らかさと当てはまりのバランスを調整する。

三角補間

cardinalTrigonometricInterpolate

template<concepts::OrderedField T>
T cardinalTrigonometricInterpolate(
    const std::vector<T>& y,
    T period,
    T xi);

動作: 等間隔の周期データ $y[0..N-1]$ (標本点 $x_k = k \cdot \text{period}/N$) を三角多項式で補間する。 Dirichlet 核アプローチ $S(x) = \frac{1}{N}\sum_k y_k\, D_N(x - x_k)$ を直接計算する (FFT 不使用、 中小規模データ向き)。 $N$ の偶奇で核を切り替える。

パラメータ:

用途: 周期信号・スペクトル法の基礎となる周期データの再構成。 標本点で厳密に元の値を再現する。

// 例: 周期 2π の正弦波を 8 点で標本化して補間
std::vector<double> y(8);
for (int k = 0; k < 8; ++k) y[k] = std::sin(2.0 * M_PI * k / 8.0);
double v = cardinalTrigonometricInterpolate<double>(y, 2.0 * M_PI, 0.3);
// 実行結果: v ≒ sin(0.3) (標本間も三角多項式で滑らかに再構成)

スプラインフィッティング (平滑化)

標本点を厳密に通す補間ではなく、 平滑化パラメータ $s$ で「当てはまり」と「滑らかさ」を調整する最小二乗フィット。 FITPACK (Dierckx) の splrep/splev に相当する。 1D・パラメトリック曲線・2D 曲面の 3 種。

関数詳細

SplineFitResult / splineFit

template<concepts::OrderedField T>
struct SplineFitResult {
    std::vector<T> knots;        // 節点ベクトル
    std::vector<T> coefficients; // B-スプライン係数
    std::size_t degree;          // スプライン次数
    T smoothing;                 // 実際の平滑化パラメータ
    T residual;                  // 残差二乗和
};

template<concepts::OrderedField T>
SplineFitResult<T> splineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& w = {},
    T s = T(-1),
    std::size_t degree = 3);

動作: データ $(x_i, y_i)$ に対し $\sum_i w_i (y_i - S(x_i))^2 + s \int (S''(x))^2 dx$ を最小化する平滑化 B-スプラインフィット。 $s = 0$ なら補間 (全点を通る)、 $s > 0$ なら平滑化 (ノイズ抑制)。 平滑化パラメータに応じて節点数を自動選択する。

パラメータ:

注意: 点数が $\text{degree}+1$ 未満、 サイズ不一致、 重み長不正では std::invalid_argument を送出。 結果は splineEval で評価する。

splineEval

template<concepts::OrderedField T>
T splineEval(const SplineFitResult<T>& fit, T x);

template<concepts::OrderedField T>
std::vector<T> splineEval(const SplineFitResult<T>& fit, std::span<const T> xs);

動作: splineFit の結果を任意点で評価する (FITPACK splev 相当)。 単点版とまとめて評価する多点版がある。

パラメータ:

// ノイズ入りデータを平滑化フィット
std::vector<double> xs = /* 昇順の x */;
std::vector<double> ys = /* ノイズ入り y */;
auto fit = splineFit<double>(xs, ys);          // s 自動、 3 次
double v = splineEval<double>(fit, 1.5);        // 平滑化曲線上の値
// fit.residual に残差二乗和、 fit.smoothing に採用された s が入る

ParametricSplineResult / parametricSplineFit

template<concepts::OrderedField T>
struct ParametricSplineResult {
    SplineFitResult<T> x_fit;   // x(t) のフィット
    SplineFitResult<T> y_fit;   // y(t) のフィット
    std::vector<T> t;            // パラメータ値
};

template<concepts::OrderedField T>
ParametricSplineResult<T> parametricSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T s = T(-1),
    std::size_t degree = 3);

動作: 平面曲線 $(x(t), y(t))$ を弧長パラメータ (累積弦長を $[0, 1]$ に正規化) でフィットする。 $x(t)$ と $y(t)$ をそれぞれ splineFit で個別にフィットする。 $x$ が単調でなくてもよい (閉曲線・ループ可)。

パラメータ:

注意: 点数 2 未満・サイズ不一致では std::invalid_argument を送出。 評価は parametricSplineEval を使う。

parametricSplineEval

template<concepts::OrderedField T>
std::pair<T, T> parametricSplineEval(
    const ParametricSplineResult<T>& fit,
    T t);

動作: パラメータ $t \in [0, 1]$ における曲線上の点 $(x(t), y(t))$ を std::pair<T, T> で返す。

パラメータ:

SurfaceSplineResult / surfaceSplineFit

template<concepts::OrderedField T>
struct SurfaceSplineResult {
    std::vector<T> knots_x;      // x 方向の節点ベクトル
    std::vector<T> knots_y;      // y 方向の節点ベクトル
    std::vector<T> coefficients; // 係数 (nx × ny の行優先フラット配列)
    std::size_t nx;             // x 方向の基底数
    std::size_t ny;             // y 方向の基底数
    std::size_t degree;         // スプライン次数
};

template<concepts::OrderedField T>
SurfaceSplineResult<T> surfaceSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& z,
    std::size_t degree = 3);

動作: 格子データ $z[i][j] = f(x_i, y_j)$ をテンソル積 B-スプライン $S(x, y) = \sum_i \sum_j c_{ij} B_i(x) B_j(y)$ でフィットする。 まず各 $x$ 行を $y$ 方向にフィットし、 続いてその係数を $x$ 方向にフィットする 2 段階分解。

パラメータ:

注意: 各方向の点数が $\text{degree}+1$ 未満、 または $z$ のサイズが $m_x m_y$ でないときは std::invalid_argument を送出。 評価は surfaceSplineEval を使う。

surfaceSplineEval

template<concepts::OrderedField T>
T surfaceSplineEval(const SurfaceSplineResult<T>& fit, T x, T y);

動作: テンソル積スプライン $S(x, y) = \sum_i \sum_j c_{ij} B_i(x) B_j(y)$ を点 $(x, y)$ で評価する。

パラメータ:

画像補間 (2D リサンプリング)

行優先の 2 次元グリッド (画像) を、 連続座標 $(x, y)$ で補間する関数群。 $x$ は列方向、 $y$ は行方向の 0 始まり座標。 境界はクランプ。

関数詳細

nearestNeighbor2D

template<concepts::OrderedField T>
T nearestNeighbor2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y);

動作: 最も近い格子点の値をそのまま返す (座標を四捨五入してインデックス化、 範囲はクランプ)。 画像の縮小やピクセルアートの拡大に使う。

パラメータ:

bicubicInterpolate

template<concepts::OrderedField T>
T bicubicInterpolate(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y);

動作: 近傍 $4 \times 4 = 16$ 格子点を Catmull-Rom カーネル ($a = -0.5$) で重み付けして滑らかに補間する。 OpenCV の INTER_CUBIC・Photoshop のバイキュービック法に相当する。 境界はインデックスをクランプ。

パラメータ:

用途: 写真の拡大・縮小で最近傍・双線形より高品質な結果が欲しいとき。

lanczosInterpolate2D

template<concepts::OrderedField T>
T lanczosInterpolate2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y,
    int a = 3);

動作: sinc を基にした高品質リサンプリング。 Lanczos-$a$ カーネル $L(x) = \mathrm{sinc}(x)\,\mathrm{sinc}(x/a)$ ($|x| < a$、 それ以外 0) を $2a \times 2a$ の窓で適用する。 境界で重み和が崩れないよう正規化する。 重みは精度を落とさないよう $T$ 型で計算する (多倍長 Float でも精度を保つ)。

パラメータ:

用途: 高品質な画像縮小・拡大。 双 3 次よりリンギングと鮮鋭度のバランスがよく、 写真リサンプリングの定番。

使用例

#include <math/interpolation/Interpolation.hpp>
#include <math/interpolation/BSpline.hpp>
#include <iostream>
#include <vector>
using namespace sangi;

int main() {
    // ---- 1) cubicSpline で点列を補間して eval する ----
    std::vector<double> xs{0.0, 1.0, 2.0, 3.0, 4.0};
    std::vector<double> ys{0.0, 1.0, 4.0, 9.0, 16.0};   // y = x^2 の標本
    auto spline = cubicSpline<double>(xs, ys);
    std::cout << "spline(2.5) = " << spline.eval(2.5) << '\n';
    // 出力: spline(2.5) ≒ 6.25 付近 (2.5^2 = 6.25)

    // ---- 2) cubic_spline_coefficients / evaluate (係数配列版) ----
    auto coeffs = cubic_spline_coefficients<double>(xs, ys);  // 自然境界
    std::cout << "eval(1.5) = "
              << cubic_spline_evaluate<double>(xs, coeffs, 1.5) << '\n';

    // ---- 3) splineFit で平滑化 (ノイズ抑制) ----
    std::vector<double> xn{0,1,2,3,4,5,6,7,8,9};
    std::vector<double> yn{0.1,0.9,2.1,2.9,4.2,4.8,6.1,6.9,8.2,8.9};  // y≒x + ノイズ
    auto fit = splineFit<double>(xn, yn);   // s 自動、 3 次
    std::cout << "fit(4.5) = " << splineEval<double>(fit, 4.5) << '\n';
    std::cout << "residual = " << fit.residual << '\n';

    // ---- 4) PCHIP で単調保存補間 ----
    std::vector<double> xm{0,1,2,3}, ym{0,0,1,1};   // 段差状の単調データ
    auto pc = pchipCoefficients<double>(xm, ym);
    std::cout << "pchip(1.5) = "
              << cubic_spline_evaluate<double>(xm, pc, 1.5) << '\n';
    // 単調保存ゆえ 0 〜 1 の範囲に収まる (オーバーシュートなし)
}

関連する数学的背景

以下の記事では、補間モジュールの基盤となる数学的概念を解説している。

• 補間とは — 標本点から間の値を再構成する基本的な考え方
• ラグランジュ補間多項式 — 多項式補間の出発点と Runge 現象
• 3 次スプライン — $C^2$ 連続な区分 3 次曲線と三重対角系
• B-スプライン — Cox-de Boor 基底と局所性
• エルミート補間 — 導関数も合わせ込む補間
• 最小二乗法 — スプラインフィッティング・回帰の基礎