エルミート曲線

この記事は, 旧ブログから移植された記事です. よって, その内容として, 旧ブログに依存した文脈が含まれている可能性があります. 予めご了承下さい.

典型的なパラメトリック曲線の一種である, エルミート曲線についてのメモ.

パラメトリック曲線

そもそもパラメトリック曲線とは, 任意のパラメータから各々の座標を陽関数形式で表現できる曲線のことをいう. このとき定義できる関数 ff がパラメータ tt1 の多項式である場合, それを多項式曲線といい, 有理式である場合, それを有理曲線という. 例えば, 直線の方程式 y=m(xa)+by = m(x-a)+b は,

{x=ty=m(ta)+b\begin{cases} x=t \\ y=m(t-a)+b \end{cases}

とパラメタライズできる. この方程式では, パラメタライズせずとも, xx に 1 つの実数を代入すれば, 必ず yy が求まる(逆も言える)ことは明らかである. 次に, 3 次曲線 y2=x3+x2y^{2} = x^{3} + x^{2} について考える. xx に 1 つの実数を代入すると, yy は,

  • x3+x2>0x^{3} + x^{2} \gt 0 のとき ±x3+x2\pm{\sqrt{x^{3} + x^{2}}},
  • x3+x2=0x^{3} + x^{2} = 0 のとき 00,
  • x3+x2<0x^{3} + x^{2} \lt 0 のとき複素数

となり, xx の値によっては, グラフが存在しない. これを

{x=t21y=t3t\begin{cases} x=t^{2}-1 \\ y=t^{3}-t \end{cases}

とパラメタライズすることで, 任意の実数 tt に対するグラフ上の一点を定めることができる2.

パラメトリック曲線との視覚的差異

エルミート曲線

エルミート曲線は, 始点 (x0,y0)(x_{0}, y_{0}) とその速度ベクトル (v0,w0)(v_{0}, w_{0}), 終点 (x1,y1)(x_{1}, y_{1}) とその速度ベクトル (v1,w1)(v_{1}, w_{1}) を与えたときに, 始点と終点の間をつなぐ 3 次パラメトリック曲線である. 3 次パラメトリック曲線と言うからには, これを一般の 3 次多項式で表現できるはずである.f(t)=At3+Bt2+Ct+D f(t) = At^{3} + Bt^{2} + Ct + D パラメトリック曲線は xxyy をそれぞれ別々に扱えるので, まず xx について考える. x(t)=At3+Bt2+Ct+D(1)x(t) = At^{3} + Bt^{2} + Ct + D \tag{1} 4 つの係数, A,B,C,DA, B, C, D が得られれば, 任意の tt に対するエルミート曲線の xx が得られるはずである. まず, 始点の座標を求める. t=0t = 0 のときを始点, t=1t = 1 のときを終点としたとき, t=0t = 0 のときに始点の座標, t=1t = 1 のときに終点の座標が得られることは明らかなので, tt00, および tt11 を代入して次の二式が得られる. x(0)=D=x0(2)x(0) = D = x_{0} \tag{2} x(1)=A+B+C+D=x1(3)x(1) = A + B + C + D = x_{1} \tag{3} ここで, ある地点 tt での曲線の傾きを求めるために, 一階微分した次の式を得る. ddtx(t)=3At2+2Bt+C(4)\displaystyle\dfrac{d}{dt}x(t)=3At^{2} + 2Bt + C \tag{4} 始点 (00) と終点 (11) をそれぞれ式 (4)(4)tt に代入し, 始点における曲線の傾きと, 終点における曲線の傾きが得られる.

dxdt(0)=C=v0(5)\displaystyle \dfrac{dx}{dt}(0) = C = v_{0} \tag{5} dxdt(1)=3A+2B+C=v1(6)\displaystyle \dfrac{dx}{dt}(1) = 3A + 2B + C = v_{1} \tag{6}

(3)(3) に式 (2)(2), (5)(5) を, 式 (6)(6) に式 (5)(5) を代入すると, 次の二式が得られる.

A+B+v0+x0=x1(7)A+B+v_{0}+x_{0} = x_{1} \tag{7} 3A+2B+v0=v1(8)3A + 2B + v_{0} = v_{1} \tag{8}

(7)(7), (8)(8) の連立方程式を解くと,

A=2x02x1+v0+v1(9)A = 2x_{0} - 2x_{1} + v_{0} + v_{1} \tag{9} B=3x0+3x12v0v1(10)B = -3x_{0}+3x_{1} - 2v_{0} - v_{1} \tag{10}

となる. 式 (2),(5),(9),(10)(2), (5), (9), (10) を式 (1)(1) に代入すると x(t)=(2x02x1+v0+v1)t3+(3x0+3x12v0v1)t2+v0t+x0(11) x(t) = (2x_{0} - 2x_{1} + v_{0} + v_{1})t^{3} + (-3x_{0} + 3x_{1} - 2v_{0} - v_{1})t^{2} + v_{0}t + x_{0} \tag{11} であるから, 式 (11)(11) から tt に対するエルミート曲線の xx が取れることがわかった. yy についても, xxyy に, vvww にすると, 同様にして得られるから

{x(t)=(2x02x1+v0+v1)t3+(3x0+3x12v0v1)t2+v0t+x0y(t)=(2y02y1+w0+w1)t3+(3y0+3y12w0w1)t2+w0t+y0 \begin{aligned} \begin{cases} x(t) &= (2x_{0} - 2x_{1} + v_{0} + v_{1})t^{3} + (-3x_{0} + 3x_{1} - 2v_{0} - v_{1})t^{2} + v_{0}t + x_{0} \\ y(t) &= (2y_{0} - 2y_{1} + w_{0} + w_{1})t^{3} + (-3y_{0} + 3y_{1} - 2w_{0} - w_{1})t^{2} + w_{0}t + y_{0} \end{cases} \end{aligned}

がいえる.

さらに, x(t)x(t) から x0,x1,v0,v1x_{0}, x_{1}, v_{0}, v_{1} (y0,y1,w0,w1y_{0}, y_{1}, w_{0}, w_{1} についても同様) の係数にそれぞれ着目して, 次のように定義する.

  • x0x_{0} の係数に着目し, 2t33t2+1=(2t+1)(1t)22t^3-3t^2+1=(2t+1)(1-t)^2, これを H03(t)H^3_{0}(t) とする
  • v0v_{0} の係数に着目し, t32t2+t=t(1t)2t^3-2t^2+t=t(1-t)^2, これを H13(t)H^3_{1}(t) とする
  • v1v_{1} の係数に着目し, t3t2=t2(1t)t^3-t^2=-t^2(1-t), これを H23(t)H^3_{2}(t) とする
  • x1x_{1} の係数に着目し, 3t22t3=t2(32t)3t^2-2t^3=t^2(3-2t), これを H33(t)H^3_{3}(t) とする

すると, 先ほど導出したx(t),y(t)x(t), y(t) の式は次のように定義できる.

エルミート曲線

{x(t)=x0H03(t)+v0H13+v1H23(t)+x1H33(t)y(t)=y0H03(t)+w0H13+w1H23(t)+y1H33(t)\begin{aligned} \begin{cases} x(t) &= x_{0}H^3_{0}(t) + v_{0}H^3_{1} + v_{1}H^3_{2}(t) + x_{1}H^3_{3}(t) \\ y(t) &= y_{0}H^3_{0}(t) + w_{0}H^3_{1} + w_{1}H^3_{2}(t) + y_{1}H^3_{3}(t) \end{cases} \end{aligned}

これは, 3 次エルミート多項式によるエルミート曲線の定義と同義である.

実際に描く

始点 (1,0)(1, 0) での速度ベクトルを (0,1)(0, 1), 終点 (1,0)(-1, 0) での速度ベクトルを (0,1)(0, -1) として, 上記に導いたエルミート曲線に従い, 点を打ってみた3.

{-# OPTIONS_GHC -Wall #-}

module Main where

import Graphics.Rendering.OpenGL 
import Graphics.UI.GLUT 
import Control.Monad (forM_)
import Control.Arrow ((&&&))

type Float' = GLfloat

linspaceWithDensity :: Float' -> Int -> Int -> [Float']
linspaceWithDensity density bt et = let distance = round (realToFrac (abs et + abs bt) / density) in 
    take distance $ iterate (+density) (realToFrac bt :: Float')

-- | The function that generates a coordinate list of Hermite curve according to
--
-- the start point, the velocity vector of the start point, the end point, the velocity vector of the end point,
-- the density of the points and the range of @t@ (@bt <= t <= et@).
hermite :: (Float', Float') -> (Float', Float') -> (Float', Float') -> (Float', Float') -> Float' -> Int -> Int -> [(Float', Float')]
hermite st stVec ed edVec = ((map (herX &&& herY).).).linspaceWithDensity
    where
        h30 t = (2 * t + 1) * (1 - t)^2
        h31 t = t * (1 - t)^2
        h32 t = -t^2 * (1 - t)
        h33 t = t^2 * (3 - 2 * t)
        hermite' t v1 v2 v3 v4 = h30 t * v1 + h31 t * v2 + h32 t * v3 + h33 t * v4
        herX t = hermite' t (fst st) (fst stVec) (fst edVec) (fst ed)
        herY t = hermite' t (snd st) (snd stVec) (snd edVec) (snd ed)

resize :: Size -> IO ()
resize s@(Size w h) = do
    viewport $= (Position 0 0, s)
    loadIdentity
    ortho (-w') w' (-h') h' (-1.0) 1.0
    where
        w' = realToFrac w / 200.0
        h' = realToFrac h / 200.0

disp :: IO ()
disp = do
    clear [ColorBuffer]
    color (Color3 0 0 0 :: Color3 GLdouble)
    pointSize $= 1.0
    renderPrimitive Points $ forM_ hcurve (vertex.(uncurry Vertex2))
    flush
    where
        hcurve = hermite (1, 0) (0, 1) (-1, 0) (0, -1) 0.001 (-2) 2

main :: IO ()
main = do
    (progname, _) <- getArgsAndInitialize
    initialDisplayMode $= [RGBAMode]
    _ <- createWindow progname
    clearColor $= Color4 1 1 1 1
    windowTitle $= "Hermite curve"
    displayCallback $= disp
    reshapeCallback $= Just resize
    mainLoop

次の曲線を得た.

エルミート曲線

  1. tt は time からくる.↩︎

  2. グラフ描画のソースコード.↩︎

  3. 浮動小数点数の絶対誤差が気になる演算だが, そのような観点において上記コードでは手抜きをした.↩︎