【S02】3次スプライン補間(変数の正規化ver.)

　
この記事はアドカレに参加しています。

スプライン補間

　3次スプライン補間とはn番目とn+1番目の点の間をn-1番目からn+2番目の4つの点で三次関数として補間するものでした。

　実は3次スプライン補間にも様々な種類が存在し、使用用途によって端の加速度を変えたりすることがあります。マイナー処理です。

変数を正規化したver.のスプライン補間もそんなマイナーの1つですが、行列が綺麗なので僕は好きです。ベジェ曲線との互換性を持とうとしたら、こちらの方が便利だと思います。ということで、変数を正規化したver.のスプライン補間の紹介です。

三次関数の式から係数を求める(変数の正規化ver.)

　まずは点の数をm+1個として、$p_{0}$～$p_{m}$という点があったとします。このとき、$p_{n}$と$p_{n+1}$を補間する三次関数$S_{n}(t)$は以下のように表すことができます。

\begin{align}
y&=S_{n}(t)　　(0 \leqq t \leqq 1)\\
y&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n}
\end{align}

$t$というのは三次関数の引数で、$t=0$のとき$S_{n}(0)=p_{n}$であり、$t=1$のとき$S_{n}(1)=p_{n+1}$となります。

係数d

係数$d_{n}$の求め方

$t=0$のときを考えます。

\begin{align}
y&=S_{n}(t) \\
p_{n}&=S_{n}(0)
\end{align}

という関係を用いて、

\begin{align}
p_{n}&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
p_{n}&=a_{n} \times 0^3+b_{n} \times 0^2+c_{n} \times 0+d_{n} \\
p_{n}&=d_{n}
\end{align}

よって、$d_{n}=p_{n}$となります。

係数c

係数$c_{n}$の求め方

関数$S_{n}(t)$の終点($p_{n+1}$)と$S_{n+1}(t)$の始点($p_{n+1}$)の値は同じです。

$t_{n}=1$、$t_{n+1}=0$として

\begin{align}
S_{n}(t_{n})&=S_{n+1}(t_{n+1}) \\
a_{n}t_{n}^3+b_{n}t_{n}^2+c_{n}t_{n}+d_{n}&=a_{n+1}t_{n+1}^3+b_{n+1}t_{n+1}^2+c_{n+1}t_{n+1}+d_{n+1} \\
a_{n} \times 1+b_{n} \times  1+c_{n} \times 1+d_{n}&=a_{n+1} \times 0+b_{n+1} \times 0+c_{n+1} \times 0+d_{n+1} \\
a_{n}+b_{n}+c_{n}+d_{n}&=d_{n+1} \\
c_{n}&=d_{n+1}-a_{n}-b_{n}-d_{n}
\end{align}

よって、$c_{n}=d_{n+1}-a_{n}-b_{n}-d_{n}$となります。

係数a

係数$a_{n}$の求め方

$S_{n}(t)$の二次微分$\frac{d^2S_{n}(t)}{dt^2}$を考えます。

\begin{align}
S_{n}(t)&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
\frac{dS_{n}(t)}{dt}&=3a_{n}t^2+2b_{n}t+c_{n} \\
\frac{d^2S_{n}(t)}{dt^2}&=6a_{n}t+2b_{n}
\end{align}

$S_{n}(t)$の始点と$S_{n+1}(t)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、$t_{n}=1$、$t_{n+1}=0$として、

\begin{align}
\frac{d^2S_{n}(t_{n})}{dt^2}&=\frac{d^2S_{n+1}(t_{n+1})}{dt^2} \\
6a_{n}t_{n}+2b_{n}&=6a_{n+1}t_{n+1}+2b_{n+1} \\
6a_{n} \times 1+2b_{n}&=6a_{n+1} \times 0+2b_{n+1} \\
6a_{n}+2b_{n}&=2b_{n+1} \\
6a_{n}&=2b_{n+1}-2b_{n} \\
a_{n}&=\frac{b_{n+1}-b_{n}}{3}
\end{align}

よって、$a_{n}=\frac{b_{n+1}-b_{n}}{3}$となります。

係数b

係数$b_{n}$の求め方

$S_{n}(t)$の一次微分$\frac{dS_{n}(t)}{dt}$を考えます。

\begin{align}
S_{n}(t)&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
\frac{dS_{n}(t)}{dt}&=3a_{n}t^2+2b_{n}t+c_{n}
\end{align}

$S_{n}(t)$の始点と$S_{n+1}(t)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、$t_{n}=1$、$t_{n+1}=0$として、

\begin{align}
\frac{dS_{n}(t_{n})}{dt}&=\frac{dS_{n+1}(t_{n+1})}{dt} \\
3a_{n}t_{n}^2+2b_{n}t_{n}+c_{n}&=3a_{n+1}t_{n+1}^2+2b_{n+1}t_{n+1}+c_{n+1} \\
3a_{n} \times 1+2b_{n}\times 1+c_{n}&=3a_{n+1}\times 0+2b_{n+1}\times 0+c_{n+1} \\
3a_{n}+2b_{n}+c_{n}&=c_{n+1} \\
\end{align}

ここで、$c_{n}=d_{n+1}-a_{n}-b_{n}-d_{n}$、$c_{n+1}=d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}$となるので、

\begin{align}
3a_{n}+2b_{n}+c_{n}&=c_{n+1} \\
3a_{n}+2b_{n}+(d_{n+1}-a_{n}-b_{n}-d_{n})&=(d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}) \\
2a_{n}+b_{n}+d_{n+1}-d_{n}&=d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}
\end{align}

項を整理します。

2a_{n}+a_{n+1}+b_{n}+b_{n+1}=d_{n}-2d_{n+1}+d_{n+2}

更に$a_{n}=\frac{b_{n+1}-b_{n}}{3}$、$a_{n+1}=\frac{b_{n+2}-b_{n+1}}{3}$として、

\begin{align}
2a_{n}+a_{n+1}+b_{n}+b_{n+1}&=d_{n}-2d_{n+1}+d_{n+2} \\
2(\frac{b_{n+1}-b_{n}}{3})+(\frac{b_{n+2}-b_{n+1}}{3})+b_{n}+b_{n+1}&=d_{n}-2d_{n+1}+d_{n+2}
\end{align}

分母が邪魔なので、両辺を3倍します。

\begin{align}
(2b_{n+1}-2b_{n})+(b_{n+2}-b_{n+1})+3b_{n}+3b_{n+1}&=3d_{n}-6d_{n+1}+3d_{n+2} \\
b_{n}+4b_{n+1}+b_{n+2}&=3d_{n}-6d_{n+1}+3d_{n+2}
\end{align}

さて、$b_{n}+4b_{n+1}+b_{n+2}=3d_{n+2}-6d_{n+1}+3d_{n}$という関係式は、あるものを彷彿とさせます。見えてきませんか？視えてきますよね？笑

$b_{n}$は二次微分なので、加速度と捉えることができます。よって、$p_{0}$と$p_{m}$の加速度を0とします。($b_{0}=b_{m}=0$)

$b_{n}$について以下の連立方程式が成り立ちます。

\begin{align}
b_{0}&=0 \\
b_{0}+4b_{1}+b_{2}&=3d_{0}-6d_{1}+3d_{2} \\
b_{1}+4b_{2}+b_{3}&=3d_{1}-6d_{2}+3d_{3} \\
b_{2}+4b_{3}+b_{4}&=3d_{2}-6d_{3}+3d_{4} \\
\vdots \\
b_{n-1}+4b_{n}+b_{n+1}&=3d_{n-1}-6d_{n}+3d_{n+1} \\
\vdots \\
b_{m-3}+4b_{m-2}+b_{m-1}&=3d_{m-3}-6d_{m-2}+3d_{m-1} \\
b_{m-2}+4b_{m-1}+b_{m}&=3d_{m-2}-6d_{m-1}+3d_{m} \\
b_{m}&=0 \\
\end{align}

はい、視えました。

\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
1 & 4 & 1 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 1 & 4 & 1 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 1 & 4 & 1 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\0 & \cdots & \cdots & 0 & 1 & 4 & 1 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 1 & 4 & 1 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}

\begin{bmatrix}
b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \\
\vdots \\ b_{m-2} \\ b_{m-1} \\ b_{m} 
\end{bmatrix}
=
\begin{bmatrix}
0 \\
3d_{0}-6d_{1}+3d_{2} \\
3d_{1}-6d_{2}+3d_{3} \\
3d_{2}-6d_{3}+3d_{4} \\
\vdots \\
3d_{m-3}-6d_{m-2}+3d_{m-1} \\
3d_{m-2}-6d_{m-1}+3d_{m} \\
0 
\end{bmatrix}

係数$b_{n}$はこの行列式を解くことで求めることができます。

係数の式一覧

それぞれの係数を求める式をまとめると、以下のようになります。

\begin{align}
d_{n}&=p_{n} \\
c_{n}&=d_{n+1}-a_{n}-b_{n}-d_{n} \\
b_{n-1}+4b_{n}+b_{n+1}&=3d_{n-1}-6d_{n}+3d_{n+1} \\
a_{n}&=\frac{b_{n+1}-b_{n}}{3}
\end{align}

プログラム例

　luaで申し訳ないですが、プログラム例です。luaの配列は1から始まるので注意が必要ですね。

--------------------------------------------------
--a,b,c,dを求める関数
--------------------------------------------------
local function spline_M(y)

local max_p=#y

local a={}
local b={}
local c={}
local d=y

b[1]=0
b[max_p]=0
for n=2,max_p-1 do
 b[n]=3*(d[n-1]-2*d[n]+d[n+1])
end

--行列の下三角を0にする
local w={0}
for i=2,max_p do
 local tmp=4-w[i-1]
 b[i]=(b[i]-b[i-1])/tmp
 w[i]=1/tmp
end

--行列の上三角を0にして、単位行列に
local i=max_p-1
while i>1 do
 b[i]=b[i]-b[i+1]*w[i]
 i=i-1
end

for n=1,max_p-1 do
 a[n]=(b[n+1]-b[n])/3
 c[n]=d[n+1]-a[n]-b[n]-d[n]
end

local sol={}
sol[1]=a
sol[2]=b
sol[3]=c
sol[4]=d

return sol
end
--------------------------------------------------
--------------------------------------------------

参考文献

・3次スプライン補間で軌跡生成：３次多項式と補間条件
・3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その１）
・3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その２）
・3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その１）
・3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その２）
・曲線描画
・C++を用いてスプライン補間を行う
・Thomas法による数値解析
・第3-2回　TDMA法（Thomas’algorithm）[python]
・３次スプライン補間

むすび

　特殊な3次スプライン補間の紹介でした。