パラメトリックスプライン曲線の計算方法

icon 項目のみ表示/展開表示の切り替え

概要

3次スプライン曲線の作成方法(パラメトリック版) xの関数となる通常のスプラインで(ノンパラメトリックスプライン曲線)は、1区間のxに対するyの値は一つしかないため、たとえば下図のx4~x5の区間のように、同じxの値に対してそれぞれ違うyを持つ曲線を描くことはできません。また、3次元空間にスプライン曲線を描画する場.は、新たな手法が必要となります。
たとえばxも別の変数から算出してしまえば、曲線の自由度が増します。このような手法をパラメトリック曲線といいます。たとえば円を角度θでパラメトリック曲線化すると、下式のようになります。
x=r \cos \theta
y=r \sin \theta

スプライン曲線も同様にx,yともに共通の変数tによって算出する

ノンパラメトリックの場.は、y=s(x)の様に表現していましたが、パラメトリックの場.はx=s(t),y=s(t)の様な式になります。この場.、tの決め方ですが、単純に0,1,2,・・・Nのような決め方が計算が楽になります。言い方を変えると、先にt=0,1,2,・・・Nと決めておき、既知点tとxに対して補間式の算定、既知点tとyに対して補間式の算定をし、x,yをそれぞれ求めて作図すればパラメトリック曲線が作図できるというわけです。補間式はノンパラメトリックy=s(x)の場.のxをtに、さらにxを求めるときはyをxに置きかえるだけです。ただしtが0,1,2,・・・Nと1ずつ増えるため、式の簡略化が可能となります。この場.、z=s(t)とすれば三次元空間にスプライン曲線を描画することができます。

既知点

下図のように6点の既知点からなる3次スプライン曲線を例に算出方法を説明する。 既知点をxn,ynとすると座標一覧は下表のとおりである。
tXY
0 0 2041.68 1575.59
1 1 2830.84 2298.99
2 2 3685.07 1900.71
3 3 4191.15 1353.67
N-1 4 4889.14 1648.74
N 5 5403.39 2536.84

ノンパラメトリック3次スプライン曲線による補間式の簡略化

3次スプライン曲線による補間式は、ノンパラメトリックスプライン曲線で算出している。ここでは、tが1置きに増加するという性質より式を簡略化する。まず、単純にノンパラメトリックスプライン曲線の式でxと表現しているものはtに置き換える。
スプライン曲線の二次導関数を算出式は、下式のとおりとなる。(算出過程は、ノンパラメトリックスプライン曲線を参照)
式を簡素化するために下式とおく。
h_n=t_{n+1}-t_n=1
h_{n+1}-h_n=t_{n+2}-t_{n+1}+t_n=x_{n+2}-x_n=2
\displaystyle k_n=6 \left( \frac{y_{n+1}-y_n}{1}-\frac{y_n-y_{n-1}}{1} \right)=6(y_{n+1}-2y_n+y_{n-1})
より一般的な表記をすると0~N-3について下記式の連立方程式が成り立つ。
n=0の時2(h_1+h_0)S_1''+h_1S_2''=k_1
n=1~N-4h_nS_n''+2(h_{n+1}+h_n)S_{n+1}''+h_{n+1}S_{n+2}''=k_{n+1}
n=N-3h_{N-3}S_{N-3}''+2(h_{n-2}+h_n)S_{n-3}''=k_{N-3}
連立方式を行列で示すと下式のとおりとなる。
\left[ \begin{matrix} 2(h_1+h_0) & h_1 & 0 \\ h_1 & 2(h_2+h_1) & h_2 \\ 0 & h_2 & 2(h_3+h_2) \\ \end{matrix} \right] \left[ \begin{matrix} S_1'' \\ S_2'' \\ S_3'' \\ \end{matrix} \right]=\left[ \begin{matrix} k_1 \\ k_2 \\ k_3 \\\end{matrix} \right]
\left[ \begin{matrix} 4 & 1& 0 \\ 1 &4 & 1 \\ 0 & 1 & 4 \\ \end{matrix} \right] \left[ \begin{matrix} S_1'' \\ S_2'' \\ S_3'' \\ \end{matrix} \right]=\left[ \begin{matrix} k_1 \\ k_2 \\ k_3 \\\end{matrix} \right]
xの補間式の算定
  t      X      hn      kn
0    0  2041.68
1    1  2830.84  1    390.42
2    2  3685.07  1   -2088.9
3    3  4191.15  1   1151.46
N-1  4  4889.14  1  -1102.44
N    5  5403.39  1
  ガウスの前進消去法を用いて上三角行列を作成する。
\left[ \begin{matrix} a11& a12& a13 \\a21& a22& a23 \\a31& a32& a33 \\ \end{matrix} \right]
上記の行列2列目以降に下式を適用する。
ui1=ai1/a11
aji = aji - ui1 * a1j
bi=bi-ui1*b1
\left[ \begin{matrix} a11& a12& a13\\ a21-u11*a11& a22-u11*a12& a23-u11*a13\\ a31-u31*a11& a32-u31*a12& a33-u31*a13 \\ \end{matrix} \right]
式が複雑なので適用結果をaij'と置く
\left[ \begin{matrix} a11'& a12'& a13'\\ 0& a22'& a23'\\ 0 &a32'& a33'\\ \end{matrix} \right]
更に3行目2列目以降の行に下式を適用する
ui2'=ai2'/a22'
aij'=aji'-ui2*a2j' bi=bi-ui1*bi
\left[ \begin{matrix} a11' &a12' &a13'\\ 0 &a22'& a23'\\ 0& a32'=a32'-ui2'*a22'& a33'=a33'-ui2'*a23'\\ \end{matrix} \right]
適用結果を整理すると下式となる
\left[ \begin{matrix} a11' & a12' & a13'  \\0&  a22'&  a23' \\ 0 & 0&  a33'=a33'-ui2'*a23' \\ \end{matrix} \right]
以下、行がなくなるまで、順番に行と列を移動しながら適用していく。
上三角行列ができるので、左上から右下方向の対角線上の値を掛けていくと行列の値が算出できる。
a11'*a22'*a33'=a33'-ui2'*a23'
前進消去で行列式を上三角行列に変換する
  1 2 3   kn
1 4 1 0  390.42
2 1 4 1 -2088.9
3 0 1 4 1151.46
ui1=ai1/a11 aji = aji - ui1 * a1j bi=bi-ui1*b1
   1     2  3       ui1
1  4     1  0     390.42
2  0  3.75  1  -2186.505  0.25
3  0     1  4    1151.46     0
aij'=aji'-ui2*a2j' bi=bi-ui1*bi
  1    2           3    kn       ui2'=ai2'/a22'
1 4    1           0               390.42 
2 0 3.75           1            -2186.505
3 0    0 3.733333333 1734.528 0.266666667
後退代入
S_3''= 1734.528 / 3.733333333 = 464.6057143
\displaystyle S_2''=\frac{ ( -2186.505 + -1 * 464.6057143 )}{ 3.75}  = -706.9629
\displaystyle
S1"~S3"の値より an,bn,cn,dnを算定する。
\displaystyle a_n=\frac{ S_{n+1}''(x)-S_n''(x))}{6(x_{n+1}-x_n)}
\displaystyle b_n=\frac{ S_n''(x))}{2}
\displaystyle c_n=\frac{(y_{n+1}-y_n)}{(x_{n+1}-x_n)}-\frac{1}{6}(x_{n+1}-x_n)(2S_n''(x_n)+S_{n+1}''(x_n))
d_n=y_n
  t      x          Sn"          a           b           c          d
0  0  2041.68             0  45.72428571          0  743.4357143  2041.68
1  1  2830.84   274.3457143  -163.551429 137.172857  880.6085714  2830.84
2  2  3685.07  -706.9628571  195.2614286 -353.48143        664.3  3685.07
3  3  4191.15   464.6057143  -77.4342857 232.302857  543.1214286  4191.15
4  4  4889.14             0            0          0       514.25  4889.14
5  5  5403.39
xを変化させながら、an,bn,cn,dnより補間値yを算出する。
y_n=S_n=a_n(x-x_n)^3+b_n(x-x_n)^2+c_n(x-x_n)+d_n
各区間のxを4分割して補間すると下表のとおりとなる。
  x      x      Sn"           a            b          c            d
0  0  2041.68               45.72428571          0  743.4357143  2041.68
                         0  45.72428571          0  743.4357143  2041.68
                         0  45.72428571          0  743.4357143  2041.68
                         0  45.72428571          0  743.4357143  2041.68
1  1  2830.84               -163.551429 137.172857  880.6085714  2830.84
               274.3457143  -163.551429 137.172857  880.6085714  2830.84
               274.3457143  -163.551429 137.172857  880.6085714  2830.84
               274.3457143  -163.551429 137.172857  880.6085714  2830.84
2  2  3685.07               195.2614286 -353.48143        664.3  3685.07
              -706.9628571  195.2614286 -353.48143        664.3  3685.07
              -706.9628571  195.2614286 -353.48143        664.3  3685.07
              -706.9628571  195.2614286 -353.48143        664.3  3685.07
3  3  4191.15               -77.4342857 232.302857  543.1214286  4191.15
               464.6057143  -77.4342857 232.302857  543.1214286  4191.15
               464.6057143  -77.4342857 232.302857  543.1214286  4191.15
               464.6057143  -77.4342857 232.302857  543.1214286  4191.15
4  4  4889.14                         0          0       514.25  4889.14
                         0            0          0       514.25  4889.14
                         0            0          0       514.25  4889.14
                         0            0          0       514.25  4889.14
5  5  4815.92  
t-tn    t    t-tn     xn
0        0     0      2041.68
      0.25  0.25  2228.253371
       0.5   0.5  2419.113393
      0.75  0.75  2618.546719
1  1     1     0      2830.84
      1.25  0.25  3057.009955
       1.5   0.5  3284.993571
      1.75  0.75  3499.457902
2  1     2     0      3685.07
      2.25  0.25  3832.103371
       2.5   0.5  3953.257321
      2.75  0.75  4066.837612
3  1     3     0      4191.15
      3.25  0.25  4340.239375
       3.5   0.5  4511.107143
      3.75  0.75  4696.493839
4  1     4     0      4889.14
      4.25  0.25    5017.7025
       4.5   0.5     5146.265
      4.75  0.75    5274.8275
5  1     5     0      5403.39
yの補間式の算定
  t       Y      hn   kn
0    0  1575.59
1    1  2298.99  1 -6730.08
2    2  1900.71  1  -892.56
3    3  1353.67  1  5052.66
N-1  4  1648.74  1  3558.18
N    5  2536.84  1
ガウスの消去法
前進消去で行列式を上三角行列に変換する
  1 2 3     kn
1 4 1 0 -6730.08
2 1 4 1  -892.56
3 0 1 4  5052.66
ui1=ai1/a11
aji = aji - ui1 * a1j bi=bi-ui1*b1
  1    2 3    ui1
1 4    1 0 -6730.08 
2 0 3.75 1   789.96 0.25
3 0    1 4  5052.66    0 
aij'=aji'-ui2*a2j' bi=bi-ui1*bi
  1    2           3 kn       ui2'=ai2'/a22'
1 4    1           0 -6730.08 
2 0 3.75           1   789.96
3 0    0 3.733333333 4842.004 0.266666667
後退代入
S_3"= 4842.004 / 3.733333333 = 1296.965357
S_2"= ( 789.96 -1 * 1296.965357 )/3.75 = -135.2014
S_1"=( -6730.08 + 0 * 1296.965357 -1 * -135.2014286 )/4 = -1648.72
S1"~S3"の値より an,bn,cn,dnを算定する。
  t    y         Sn"           a          b             c        d
 0 0 1575.59            0 -274.786607          0  998.1866071 1575.59
 1 1 2298.99 -1648.719643 252.2530357 -824.35982  173.8267857 2298.99
 2 2 1900.71 -135.2014286 238.6944643 -67.600714   -718.13375 1900.71
 3 3 1353.67  1296.965357 -216.160893 648.482679 -137.2517857 1353.67
 4 4 1648.74            0           0          0        888.1 1648.74
 5 5 2536.84
xを変化させながら、an,bn,cn,dnより補間値yを算出する。 各区間のxを4分割して補間すると下表のとおりとなる。
  y        y          Sn"           a         b             c          d
0  0  1575.59                  -274.786607          0   998.1866071  1575.59
                            0  -274.786607          0   998.1866071  1575.59
                            0  -274.786607          0   998.1866071  1575.59
                            0  -274.786607          0   998.1866071  1575.59
1  1  2298.99                  252.2530357 -824.35982   173.8267857  2298.99
                 -1648.719643  252.2530357 -824.35982   173.8267857  2298.99
                 -1648.719643  252.2530357 -824.35982   173.8267857  2298.99
                 -1648.719643  252.2530357 -824.35982   173.8267857  2298.99
2  2  1900.71                  238.6944643 -67.600714    -718.13375  1900.71
                 -135.2014286  238.6944643 -67.600714    -718.13375  1900.71
                 -135.2014286  238.6944643 -67.600714    -718.13375  1900.71
                 -135.2014286  238.6944643 -67.600714    -718.13375  1900.71
3  3  1353.67                  -216.160893 648.482679  -137.2517857  1353.67
                  1296.965357  -216.160893 648.482679  -137.2517857  1353.67
                  1296.965357  -216.160893 648.482679  -137.2517857  1353.67
                  1296.965357  -216.160893 648.482679  -137.2517857  1353.67
4  4  1648.74                            0          0         888.1  1648.74
                            0            0          0         888.1  1648.74
                            0            0          0         888.1  1648.74
                            0            0          0         888.1  1648.74
5  5  2536.84
t-tn  t  t-tn  yn  
0        0     0      1575.59
      0.25  0.25  1820.843111
       0.5   0.5  2040.334978
      0.75  0.75  2208.304355
1  1     1     0      2298.99
      1.25  0.25  2294.865661
       1.5   0.5  2211.345067
      1.75  0.75  2072.076939
2  1     2     0      1900.71
      2.25  0.25  1720.681119
       2.5   0.5  1554.579754
      2.75  0.75  1424.783513
3  1     3     0      1353.67
      3.25  0.25  1356.509707
       3.5   0.5  1420.144665
      3.75  0.75  1524.309791
4  1     4     0      1648.74
      4.25  0.25     1870.765
       4.5   0.5      2092.79
      4.75  0.75     2314.815
5  1     5     0      2536.84
PDF版