bate's blog

調べたこと実装したことなどを取りとめもなく書きます。

PythonとSciPyで計算した

Pythonを使いこなせていないが、それでも何とか計算だけさせた。
通過点
 p_{0} = ( -2, 0 ), \, p_{1} = ( -1, 1 ), \, p_{2} = (0,0),\, p_{3} = (1,-1)
求めるコントロールポイント
 q_{i}\hspace{3pt}(i=0,1,2,3)
ノットベクトル
 T= [ \,-3,-2,-1,0,1,2,3,4 \, ]
パラメータ
 t_{i} = \begin{cases} 0 & (i=0) \\ t_{i}+\frac{ | p_{i} - p_{i-1} | }{d} & (i=1,\dots,np-2) \\ 1 & (i=np-1) \end{cases}
 d = \sum_{1}^{np-1} \,|p_{i}-p_{i-1}|
 \hspace{5pt} = 3 \sqrt{ 2 }
上記で、パラメータを計算すると
 t_{0} = 0,\, t_{1} = 0+\frac{ \sqrt{2} }{ 3\sqrt{2} } = \frac{1}{3}, \, t_{2} = \frac{1}{3} + \frac{ \sqrt{2}}{3\sqrt{2}}=\frac{2}{3}, \, t_{3} = 1
de Boor Cox漸化式で基底関数を求める
 N_{0,4}(t) = \frac{1}{6}(1-t)^3, \hspace{5pt} N_{1,4}(t) = \frac{1}{2}t^3-t^2+\frac{2}{3}, \hspace{5pt} N_{2,4}(t) = -\frac{1}{2}t^3+\frac{1}{2}t^2+\frac{1}{2}t+\frac{1}{6}, \hspace{5pt} N_{3,4}(t)=\frac{1}{6}t^3
 t_{0}=0の時
 N_{0,4}(t_{0})=\frac{1}{6}, \hspace{5pt} N_{1,4}(t_{0})=\frac{2}{3},\hspace{5pt} N_{2,4}(t_{0})=\frac{1}{6}, \hspace{5pt} N_{3,4}(t_{0})=0
 t_{1}=\frac{1}{3}の時
 N_{0,4}(t_{1})=\frac{4}{81},\hspace{5pt} N_{1,4}(t_{1})=\frac{31}{54},\hspace{5pt} N_{2,4}(t_{1})=\frac{10}{27},\hspace{5pt} N_{3,4}(t_{1})=\frac{1}{162}
 t_{2} = \frac{2}{3}の時
 N_{0,4}(t_{2})=\frac{1}{162},\hspace{5pt} N_{1,4}(t_{2})=\frac{10}{27},\hspace{5pt} N_{2,4}(t_{2})=\frac{31}{54},\hspace{5pt} N_{3,4}(t_{2})=\frac{4}{81}
 t_{3}=1の時
 N_{0,4}(t_{3})=0,\hspace{5pt} N_{1,4}(t_{3})=\frac{1}{6},\hspace{5pt} N_{2,4}(t_{3})=\frac{2}{3},\hspace{5pt} N_{3,4}(t_{3})=\frac{1}{6}

 \begin{pmatrix} p_{0} \\ p_{1} \\ p_{2} \\ p_{3} \end{pmatrix} = \begin{pmatrix} N_{0,4}(t_{0}) & N_{1,4}(t_{0}) & N_{2,4}(t_{0}) & N_{3,4}(t_{0}) \\ N_{0,4}(t_{1}) & N_{1,4}(t_{1}) & N_{2,4}(t_{1}) & N_{3,4}(t_{1}) \\ N_{0,4}(t_{2}) & N_{1,4}(t_{2}) & N_{2,4}(t_{2}) & N_{3,4}(t_{2}) \\ N_{0,4}(t_{3}) & N_{1,4}(t_{3}) & N_{2,4}(t_{3}) & N_{3,4}(t_{3}) \end{pmatrix}  \begin{pmatrix} q_{0} \\ q_{1} \\ q_{2} \\ q_{3} \end{pmatrix}

 \begin{pmatrix} q_{0} \\ q_{1} \\ q_{2} \\ q_{3} \end{pmatrix} = \begin{pmatrix} N_{0,4}(t_{0}) & N_{1,4}(t_{0}) & N_{2,4}(t_{0}) & N_{3,4}(t_{0}) \\ N_{0,4}(t_{1}) & N_{1,4}(t_{1}) & N_{2,4}(t_{1}) & N_{3,4}(t_{1}) \\ N_{0,4}(t_{2}) & N_{1,4}(t_{2}) & N_{2,4}(t_{2}) & N_{3,4}(t_{2}) \\ N_{0,4}(t_{3}) & N_{1,4}(t_{3}) & N_{2,4}(t_{3}) & N_{3,4}(t_{3}) \end{pmatrix}^{-1} \begin{pmatrix} p_{0} \\ p_{1} \\ p_{2} \\ p_{3} \end{pmatrix}

 N = \begin{pmatrix} \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 \\ \frac{4}{81} & \frac{31}{54} & \frac{20}{54} & \frac{1}{162} \\ \frac{1}{162} & \frac{20}{54} & \frac{31}{54} & \frac{4}{81} \\ 0 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{pmatrix}

>>> from scipy import *
>>> a11 = 1.0/6.0
>>> a12 = 2.0/3.0
>>> a13 = a11
>>> a14 = 0.0
>>> a1 = [ a11, a12, a13, a14 ]
>>> a1
[0.16666666666666666, 0.66666666666666663, 0.16666666666666666, 0.0]
>>> a21 = 4.0/81.0
>>> a22 = 31.0/54.0
>>> a23 = 20.0/54.0
>>> a24 = 1.0/162.0
>>> a2 = [ a21, a22, a23, a24 ]
>>> a2
[0.049382716049382713, 0.57407407407407407, 0.37037037037037035, 0.0061728395061728392]
>>> a31 = 1.0/162.0
>>> a32 = 20.0/54.0
>>> a33 = 31.0/54.0
>>> a34 = 4.0/81.0
>>> a3 = [ a31, a32, a33, a34 ]
>>> a3
[0.0061728395061728392, 0.37037037037037035, 0.57407407407407407, 0.049382716049382713]
>>> a41 = 0.0
>>> a42 = 1.0/6.0
>>> a43 = 2.0/3.0
>>> a44 = 1.0/6.0
>>> a4 = [a41, a42, a43, a44]
>>> a4
[0.0, 0.16666666666666666, 0.66666666666666663, 0.16666666666666666]
>>> a = [ a1, a2, a3, a4 ]
>>> a
>>> A = mat( a )
>>> A
matrix([[ 0.16666667,  0.66666667,  0.16666667,  0.        ],
        [ 0.04938272,  0.57407407,  0.37037037,  0.00617284],
        [ 0.00617284,  0.37037037,  0.57407407,  0.04938272],
        [ 0.        ,  0.16666667,  0.66666667,  0.16666667]])
>>> B = A.I
>>> B
matrix([[ 12.5, -24. ,  16.5,  -4. ],
        [ -2. ,   7.5,  -6. ,   1.5],
        [  1.5,  -6. ,   7.5,  -2. ],
        [ -4. ,  16.5, -24. ,  12.5]])
>>> A*B
matrix([[  1.00000000e+00,  -2.22044605e-16,   2.22044605e-16,
           0.00000000e+00],
        [  1.04083409e-16,   1.00000000e+00,   2.77555756e-16,
          -2.77555756e-17],
        [  1.66533454e-16,   4.44089210e-16,   1.00000000e+00,
           0.00000000e+00],
        [  0.00000000e+00,   4.44089210e-16,  -4.44089210e-16,
           1.00000000e+00]])
>>> B*A
matrix([[  1.00000000e+00,   1.66533454e-15,   4.44089210e-16,
          -1.11022302e-16],
        [ -9.02056208e-17,   1.00000000e+00,  -2.22044605e-16,
           0.00000000e+00],
        [  1.59594560e-16,   7.77156117e-16,   1.00000000e+00,
           0.00000000e+00],
        [ -3.88578059e-16,  -1.77635684e-15,   0.00000000e+00,
           1.00000000e+00]])

>>> #通過点
>>> p0 = [ -2.0, 0.0 ]
>>> p1 = [ -1.0, 1.0 ]
>>> p2 = [ 0.0, 0.0 ]
>>> p3 = [ 1.0, -1.0 ]

>>> #コントロールポイントを求める
>>> p_x = [ p0[0], p1[0], p2[0], p3[0] ]
>>> p_x
[-2.0, -1.0, 0.0, 1.0]
>>> (B[0]).A * p_x
array([[-25.,  24.,   0.,  -4.]])
>>> q0_x = -25. + 24. + 0. -4.
>>> q0_x
-5.0
>>> (B[1]).A * p_x
array([[ 4. , -7.5,  0. ,  1.5]])
>>> q1_x = 4. - 7.5 + 0. + 1.5
>>> q1_x
-2.0
>>> (B[2]).A * p_x
array([[-3.,  6.,  0., -2.]])
>>> q2_x = -3. + 6. + 0. -2.
>>> q2_x
1.0
>>> (B[3]).A * p_x
array([[  8. , -16.5,   0. ,  12.5]])
>>> q3_x = 8. -16.5 + 0. + 12.5
>>> q3_x
4.0

>>> p_y = [ p0[1], p1[1], p2[1], p3[1] ]
>>> p_y
[0.0, 1.0, 0.0, -1.0]
>>> (B[0]).A * p_y
array([[  0., -24.,   0.,   4.]])
>>> q0_y = 0. -24. + 0. + 4.
>>> q0_y
-20.0
>>> q0 = [ q0_x, q0_y ]
>>> q0
[-5.0, -20.0]
>>> (B[1]).A * p_y
array([[ 0. ,  7.5,  0. , -1.5]])
>>> q1_y = 0. + 7.5 + 0. -1.5
>>> q1_y
6.0
>>> q1 = [ q1_x, q1_y ]
>>> q1
[-2.0, 6.0]
>>> (B[2]).A * p_y
array([[ 0., -6.,  0.,  2.]])
>>> q2_y = 0. -6. + 0. + 2.
>>> q2_y
-4.0
>>> q2 = [ q2_x, q2_y ]
>>> q2
[1.0, -4.0]
>>> (B[3]).A * p_y
array([[  0. ,  16.5,   0. , -12.5]])
>>> q3_y = 0. + 16.5 + 0. -12.5
>>> q3_y
4.0
>>> q3 = [ q3_x, q3_y ]
>>> q3
[4.0, 4.0]

>>> #通過点を通るか確認
>>> 1.0/6.0*q0[0] + 2.0/3.0*q1[0] + 1.0/6.0*q2[0]
-1.9999999999999998
>>> 1.0/6.0*q0[1] + 2.0/3.0*q1[1] + 1.0/6.0*q2[1]
3.3306690738754696e-016
>>> p0
[-2.0, 0.0]
>>> 4.0/81.0*q0[0] + 31.0/54.0*q1[0] + 10.0/27.0*q2[0] + 1.0/162.0*q3[0]
-1.0
>>> 4.0/81.0*q0[1] + 31.0/54.0*q1[1] + 10.0/27.0*q2[1] + 1.0/162.0*q3[1]
1.0000000000000004
>>> p1
[-1.0, 1.0]
>>> 1.0/162.0*q0[0] + 10.0/27.0*q1[0] + 31.0/54.0*q2[0] + 4.0/81.0*q3[0]
0.0
>>> 1.0/162.0*q0[1] + 10.0/27.0*q1[1] + 31.0/54.0*q2[1] + 4.0/81.0*q3[1]
0.0
>>> p2
[0.0, 0.0]
>>> 1.0/6.0*q1[0] + 2.0/3.0*q2[0] + 1.0/6.0*q3[0]
1.0
>>> 1.0/6.0*q1[1] + 2.0/3.0*q2[1] + 1.0/6.0*q3[1]
-0.99999999999999989
>>> p3
[1.0, -1.0]