3次エルミートスプライン

数値解析において、3次エルミートスプラインまたは3次エルミート補間は、各部分がエルミート形式で指定される3次多項式、つまり対応する区間の両端の値と1次微分によって指定されるスプラインである ^[1] 。

3次エルミートスプラインは、通常、与えられる引数値 ${\displaystyle x_{1},x_{2},\ldots ,x_{n}}$ で指定される数値データを補間して連続関数を得るために使用される。 データは、各 ${\displaystyle x_{k}}$ についての所望の関数値と微分から構成される必要がある（値のみが与えられている場合、微分はそれから推定する必要がある）。 エルミート補間多項式^{[訳注 1]}は各区間 ${\displaystyle (x_{k},x_{k+1})}$ に個別に適用される。 結果のスプラインは連続、かつ、一次微分連続となる。

3次多項式スプラインは他の方法でも指定でき、ベジェ形式が最も一般的である。 しかし、これら2つの方法で得られるスプラインセットは同一であり、ベジェ形式とエルミート形式の間でデータの変換も容易であるため、これらの名称はしばしば同義語のように使用される。

3次多項式スプラインは、コンピュータグラフィックスや幾何モデリング（英語: geometric modeling）において、平面または3次元空間上の指定された点を通過する曲線や運動軌跡を得るために広く用いられている。 これらのアプリケーションでは、平面または空間の各座標は、個別のパラメータ t を持つ3次スプライン関数によって個別に補間される。 3次多項式スプラインは、オイラー・ベルヌーイ梁理論（英語: Euler–Bernoulli beam theory）などの構造解析アプリケーションでも広く用いられている。

3次多項式スプラインは、死亡率分析^[2] や死亡率予測^[3] にも応用されている。

3次スプラインは、いくつかの方法で2つ以上のパラメータを持つ関数に拡張できる。 バイキュービックスプライン（バイキュービック補間）は、デジタル画像のピクセル値や地形の高度データなど、規則的な矩形グリッド上のデータを補間するためによく使用される。 バイキュービック曲面パッチ ^{[訳注 2]} は、コンピュータグラフィックスにおいて不可欠なツールである。

3次スプラインは、特にコンピュータグラフィックスではしばしばcsplineと呼ばれる。エルミートスプラインは、シャルル・エルミートにちなんで名付けられた。

単位間隔 ${\displaystyle [0,1]}$ について、 ${\displaystyle t=0}$ における始点 ${\displaystyle {\boldsymbol {p}}_{0}}$、${\displaystyle t=1}$ における終点 ${\displaystyle {\boldsymbol {p}}_{1}}$、および、 ${\displaystyle t=0}$ における開始接ベクトル ${\displaystyle {\boldsymbol {m}}_{0}}$、${\displaystyle t=1}$ における終了接ベクトル ${\displaystyle {\boldsymbol {m}}_{1}}$ が与えられると、 多項式は次のように定義される。

${\displaystyle {\boldsymbol {p}}(t)=\left(2t^{3}-3t^{2}+1\right){\boldsymbol {p}}_{0}+\left(t^{3}-2t^{2}+t\right){\boldsymbol {m}}_{0}+\left(-2t^{3}+3t^{2}\right){\boldsymbol {p}}_{1}+\left(t^{3}-t^{2}\right){\boldsymbol {m}}_{1}}$

ただし t ∈ [0, 1] である。

任意の区間 ${\displaystyle (x_{k},x_{k+1})}$ 内の ${\displaystyle x}$ での補間は、アフィン（1次）変数変換により区間を ${\displaystyle [0,1]}$ にマッピングすることで行われる。 式は次のようになる。

${\displaystyle {\boldsymbol {p}}(x)=h_{00}(t){\boldsymbol {p}}_{k}+h_{10}(t)(x_{k+1}-x_{k}){\boldsymbol {m}}_{k}+h_{01}(t){\boldsymbol {p}}_{k+1}+h_{11}(t)(x_{k+1}-x_{k}){\boldsymbol {m}}_{k+1}}$

ただし ${\displaystyle t=(x-x_{k})/(x_{k+1}-x_{k})}$ と ${\displaystyle h}$ は、以下に定義される基底関数を参照。 なお接ベクトルの長さは、単位区間の方程式に対し ${\displaystyle x_{k+1}-x_{k}}$ 倍にスケーリングされる。

上記の式は、指定された接ベクトルを持つ2点間の、一意の3次多項式の経路を示す。

${\displaystyle P,Q}$ は、与えられた境界条件を満たす2つの3次多項式とする。${\displaystyle R=Q-P}$ とすると、

${\displaystyle R(0)=Q(0)-P(0)=0}$

${\displaystyle R(1)=Q(1)-P(1)=0}$

である。 ${\displaystyle Q}$ と ${\displaystyle P}$ はいずれも3次多項式であるので、${\displaystyle R}$ は最大で3次の多項式である。 そのため、${\displaystyle R}$ は以下の式でなければならない。

${\displaystyle R(x)=ax(x-1)(x-r)}$

導関数を計算すると、

${\displaystyle R'(x)=ax(x-1)+ax(x-r)+a(x-1)(x-r)}$

さらに、以下が明らかである。

${\displaystyle R'(0)=Q'(0)-P'(0)=0,}$

${\displaystyle R'(0)=0=ar}$

(1)

${\displaystyle R'(1)=Q'(1)-P'(1)=0,}$

${\displaystyle R'(1)=0=a(1-r)}$

(2)

式 (1) と (2) を合わせると、${\displaystyle a=0}$ が導かれ、したがって ${\displaystyle R=0}$、よって ${\displaystyle P=Q}$ となる。

単位区間上の補間多項式は次のように書ける（任意区間については、上記再スケール版を参照）。 $${\displaystyle {\boldsymbol {p}}(t)=h_{00}(t){\boldsymbol {p}}_{0}+h_{10}(t){\boldsymbol {m}}_{0}+h_{01}(t){\boldsymbol {p}}_{1}+h_{11}(t){\boldsymbol {m}}_{1}}$$ ただし ${\displaystyle h_{00}}$, ${\displaystyle h_{10}}$, ${\displaystyle h_{01}}$, ${\displaystyle h_{11}}$ はエルミート基底関数である。 これらは様々な方法で記述することができ、それぞれ異なる性質を示す。

	展開	因数分解	バーンスタイン
${\displaystyle h_{00}(t)}$	${\displaystyle 2t^{3}-3t^{2}+1}$	${\displaystyle (1+2t)(1-t)^{2}}$	${\displaystyle B_{0}(t)+B_{1}(t)}$
${\displaystyle h_{10}(t)}$	${\displaystyle t^{3}-2t^{2}+t}$	${\displaystyle t(1-t)^{2}}$	${\displaystyle {\tfrac {1}{3}}B_{1}(t)}$
${\displaystyle h_{01}(t)}$	${\displaystyle -2t^{3}+3t^{2}}$	${\displaystyle t^{2}(3-2t)}$	${\displaystyle B_{3}(t)+B_{2}(t)}$
${\displaystyle h_{11}(t)}$	${\displaystyle t^{3}-t^{2}}$	${\displaystyle t^{2}(t-1)}$	${\displaystyle -{\tfrac {1}{3}}B_{2}(t)}$

「展開」列は上記の定義で使用された表現を示している。 「因数分解」列は、${\displaystyle h_{10}}$ と ${\displaystyle h_{11}}$ が端点でゼロであることを直接示している。 さらに、${\displaystyle h_{01}}$ と ${\displaystyle h_{11}}$ は0で重複度2の零点を持ち、${\displaystyle h_{00}}$ と ${\displaystyle h_{10}}$ は1で重複度2の零点を持つため、それら端点で傾きが0になると言える。 「バーンスタイン」列は、エルミート基底関数を3次のバーンスタイン多項式に分解したものを示しており、

${\displaystyle B_{k}(t)={\binom {3}{k}}\cdot t^{k}\cdot (1-t)^{3-k}}$

である。

この関係を利用すると、3次エルミート補間を 3次ベジェ曲線の形式で、次の4つの値 ${\displaystyle {\boldsymbol {p}}_{0},{\boldsymbol {p}}_{0}+{\frac {1}{3}}{\boldsymbol {m}}_{0},{\boldsymbol {p}}_{1}-{\frac {1}{3}}{\boldsymbol {m}}_{1},{\boldsymbol {p}}_{1}}$ によって表現でき、ド・カステリョのアルゴリズムを用いてエルミート補間を行うことができる。 これは、3次ベジェ曲線 ^{[訳注 3]} において、補間曲線の両端点における接ベクトルが、中間の2つの制御点により決定されることを示している。

この多項式は次の標準形式で記述することもできる： $${\displaystyle {\boldsymbol {p}}(t)=\left(2{\boldsymbol {p}}_{0}+{\boldsymbol {m}}_{0}-2{\boldsymbol {p}}_{1}+{\boldsymbol {m}}_{1}\right)t^{3}+\left(-3{\boldsymbol {p}}_{0}+3{\boldsymbol {p}}_{1}-2{\boldsymbol {m}}_{0}-{\boldsymbol {m}}_{1}\right)t^{2}+{\boldsymbol {m}}_{0}t+{\boldsymbol {p}}_{0}}$$ ただし、制御点と接ベクトルが係数である。 この形式では係数が一定となり、一度計算して再利用できるため、t の様々な値について多項式を効率的に評価できる。

行列形式では $${\displaystyle {\boldsymbol {p}}(t)={\begin{bmatrix}t^{3}&t^{2}&t&1\end{bmatrix}}{\begin{bmatrix}2&1&-2&1\\-3&-2&3&-1\\0&1&0&0\\1&0&0&0\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {p}}_{0}\\{\boldsymbol {m}}_{0}\\{\boldsymbol {p}}_{1}\\{\boldsymbol {m}}_{1}\end{bmatrix}}}$$ となる。

3次ベジェ曲線の行列形式は $${\displaystyle {\boldsymbol {p}}(t)={\begin{bmatrix}t^{3}&t^{2}&t&1\end{bmatrix}}{\begin{bmatrix}-1&3&-3&1\\3&-6&3&0\\-3&3&0&0\\1&0&0&0\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {p_{bz}}}_{0}\\{\boldsymbol {p_{bz}}}_{1}\\{\boldsymbol {p_{bz}}}_{2}\\{\boldsymbol {p_{bz}}}_{3}\end{bmatrix}}}$$ であることから、3次エルミートスプラインと等しくなる3次ベジェ曲線の制御点は $${\displaystyle {\begin{bmatrix}{\boldsymbol {p_{bz}}}_{0}\\{\boldsymbol {p_{bz}}}_{1}\\{\boldsymbol {p_{bz}}}_{2}\\{\boldsymbol {p_{bz}}}_{3}\end{bmatrix}}={\begin{bmatrix}-1&3&-3&1\\3&-6&3&0\\-3&3&0&0\\1&0&0&0\end{bmatrix}}^{-1}{\begin{bmatrix}2&1&-2&1\\-3&-2&3&-1\\0&1&0&0\\1&0&0&0\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {p}}_{0}\\{\boldsymbol {m}}_{0}\\{\boldsymbol {p}}_{1}\\{\boldsymbol {m}}_{1}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {p}}_{0}\\{\boldsymbol {p}}_{0}+{\frac {{\boldsymbol {m}}_{0}}{3}}\\{\boldsymbol {p}}_{1}-{\frac {{\boldsymbol {m}}_{1}}{3}}\\{\boldsymbol {p}}_{1}\end{bmatrix}}}$$ である。

あるデータセット ${\displaystyle (x_{k},{\boldsymbol {p}}_{k})}$（${\displaystyle k=1,\ldots ,n}$）は、各区間に上記手順を適用することで補間できる。ただし、接ベクトルを適切な方法で選択する、つまり端点を共有する区間どうしで接線を等しくする。 補間された曲線は区分3次エルミートスプラインで構成され、${\displaystyle (x_{1},x_{n})}$ において大域的に連続して微分可能である。

接ベクトルの選択は一意ではなく、いくつかのオプションが利用可能である。

最も単純な選択肢は3点差分であり、一定間隔のデータ点を必要としない。 内部の点 ${\displaystyle k=2,\dots ,n-1}$ に対しては、

${\displaystyle {\boldsymbol {m}}_{k}={\frac {1}{2}}\left({\frac {{\boldsymbol {p}}_{k+1}-{\boldsymbol {p}}_{k}}{x_{k+1}-x_{k}}}+{\frac {{\boldsymbol {p}}_{k}-{\boldsymbol {p}}_{k-1}}{x_{k}-x_{k-1}}}\right)}$

とし、データ集合の端点では片側差分とする。

${\displaystyle {\boldsymbol {m}}_{k}=(1-c){\frac {{\boldsymbol {p}}_{k+1}-{\boldsymbol {p}}_{k-1}}{x_{k+1}-x_{k-1}}}}$

を用いて接ベクトルを計算すると、カーディナルスプライン（英語: cardinal spline）が得られる ^[4] 。なお、カーディナルスプラインは正準スプライン（英語: canonical spline）と呼ばれることもある ^[5] 。 パラメータ c はテンションパラメータであり、[0, 1] の範囲でなければならない。 ある意味では、これは接ベクトルの「長さ」と解釈できる。 c = 1 を選択するとすべての接ベクトルがゼロになり、c = 0 を選択すると一様パラメータ化（英語: uniform parameterization）Catmull-Romスプラインになる。

接ベクトルとして

${\displaystyle {\boldsymbol {m}}_{k}={\frac {{\boldsymbol {p}}_{k+1}-{\boldsymbol {p}}_{k-1}}{2}}}$

を選択することによりカーディナルスプラインの特殊な場合のひとつであるCatmull-Romスプラインが得られる。 「一様な」パラメータ間隔を前提とする。

Kochanek–Bartels スプラインは、データ ポイント ${\displaystyle {\boldsymbol {p}}_{k-1}}$, ${\displaystyle {\boldsymbol {p}}_{k}}$, ${\displaystyle {\boldsymbol {p}}_{k+1}}$ が与えられた場合に接ベクトルを選択する方法をさらに一般化したものであり、テンション、バイアス、連続性の3つのパラメータを指定できる。

上記のいずれかのタイプの3次エルミートスプラインを単調データセットの補間に使用する場合、補間された関数は必ずしも単調にはならないが、接ベクトルを調整することで単調性を維持できる。

点 ${\displaystyle {\boldsymbol {p}}_{n-1},{\boldsymbol {p}}_{n},{\boldsymbol {p}}_{n+1},{\boldsymbol {p}}_{n+2}}$ のひとつの座標成分を、 整数座標 x = n − 1, n, n + 1, n + 2 で関数 f(x) が取る値とすると、

${\displaystyle p_{n}=f(n)\quad \forall n\in \mathbb {Z} }$

と表すことができる。

さらに、端点における接ベクトルは、隣接点の中心差分 $${\displaystyle m_{n}={\frac {f(n+1)-f(n-1)}{2}}={\frac {p_{n+1}-p_{n-1}}{2}}\quad \forall n\in \mathbb {Z} .}$$ で定義されると仮定する。

実数値 x における補間値 f(x) を評価するには、まず x を整数部 n と小数部 u に分離する。

${\displaystyle x=n+u,}$

${\displaystyle n=\lfloor x\rfloor =\operatorname {floor} (x),}$

${\displaystyle u=x-n=x-\lfloor x\rfloor ,}$

${\displaystyle 0\leq u<1,}$

ここで ${\displaystyle \lfloor x\rfloor }$ は床関数を表し、x 以下の最大の整数を返すものとする。

Catmull-Romスプラインは ^[7] $${\displaystyle {\begin{aligned}f(x)=f(n+u)&={\text{CINT}}_{u}(p_{n-1},p_{n},p_{n+1},p_{n+2})\\&={\begin{bmatrix}1&u&u^{2}&u^{3}\end{bmatrix}}{\begin{bmatrix}0&1&0&0\\-{\tfrac {1}{2}}&0&{\tfrac {1}{2}}&0\\1&-{\tfrac {5}{2}}&2&-{\tfrac {1}{2}}\\-{\tfrac {1}{2}}&{\tfrac {3}{2}}&-{\tfrac {3}{2}}&{\tfrac {1}{2}}\end{bmatrix}}{\begin{bmatrix}p_{n-1}\\p_{n}\\p_{n+1}\\p_{n+2}\end{bmatrix}}\\&={\frac {1}{2}}{\begin{bmatrix}-u^{3}+2u^{2}-u\\3u^{3}-5u^{2}+2\\-3u^{3}+4u^{2}+u\\u^{3}-u^{2}\end{bmatrix}}^{\mathrm {T} }{\begin{bmatrix}p_{n-1}\\p_{n}\\p_{n+1}\\p_{n+2}\end{bmatrix}}\\&={\frac {1}{2}}{\begin{bmatrix}u{\big (}(2-u)u-1{\big )}\\u^{2}(3u-5)+2\\u{\big (}(4-3u)u+1{\big )}\\u^{2}(u-1)\end{bmatrix}}^{\mathrm {T} }{\begin{bmatrix}p_{n-1}\\p_{n}\\p_{n+1}\\p_{n+2}\end{bmatrix}}\\&={\tfrac {1}{2}}{\Big (}{\big (}u^{2}(2-u)-u{\big )}p_{n-1}+{\big (}u^{2}(3u-5)+2{\big )}p_{n}+{\big (}u^{2}(4-3u)+u{\big )}p_{n+1}+u^{2}(u-1)p_{n+2}{\Big )}\\&={\tfrac {1}{2}}{\big (}(-u^{3}+2u^{2}-u)p_{n-1}+(3u^{3}-5u^{2}+2)p_{n}+(-3u^{3}+4u^{2}+u)p_{n+1}+(u^{3}-u^{2})p_{n+2}{\big )}\\&={\tfrac {1}{2}}{\big (}(-p_{n-1}+3p_{n}-3p_{n+1}+p_{n+2})u^{3}+(2p_{n-1}-5p_{n}+4p_{n+1}-p_{n+2})u^{2}+(-p_{n-1}+p_{n+1})u+2p_{n}{\big )}\\&={\tfrac {1}{2}}{\Big (}{\big (}(-p_{n-1}+3p_{n}-3p_{n+1}+p_{n+2})u+(2p_{n-1}-5p_{n}+4p_{n+1}-p_{n+2}){\big )}u+(-p_{n-1}+p_{n+1}){\Big )}u+p_{n},\end{aligned}}}$$ ただし、${\displaystyle \mathrm {T} }$ は行列の転置を表す。 最後の式はホーナー法を用いた計算手法を示す。

この記述はトリキュービック補間（1回の最適化につき、同じ u と異なる p を使用して CINT_u を16回計算する必要がある）にも関連する。

ポータル 数学

• バイキュービック補間、2次元への一般化
• トリキュービック補間（英語: Tricubic interpolation）、3次元への一般化
• エルミート補間（英語: Hermite interpolation）
• 多変量補間（英語: Multivariate interpolation）
• スプライン補間
• 離散スプライン補間（英語: Discrete spline interpolation）

• Spline Curves, Prof. Donald H. House Clemson University
• Multi-dimensional Hermite Interpolation and Approximation, Prof. Chandrajit Bajaj, Purdue University
• Introduction to Catmull–Rom Splines, MVPs.org
• Interpolating Cardinal and Catmull–Rom splines
• Interpolation methods: linear, cosine, cubic and hermite (with C sources)
• Common Spline Equations