キャットムル-ロム スプライン曲線

最近、ごく基本的なの曲線のアルゴリズムを実装したりしてる。

キャットムル-ロム スプライン(Catmull-Rom Spline)ていうアルゴリズムは、制御点を必ず通るという特徴があるので、絵を生成するような用途だととても使いやすい。しかも実装が簡単。 最近、これを実装する機会があったので、仕組みを説明してみようと思う。

スプライン曲線

ゲームのアルゴリズムの本をめくり、曲線について書かれているページを開くと、大抵は、次数がいくつであっても対応できる級数みたいなものの式が載っていて、難しそうである。

これはどうしてかというと、数学をもってすれば、めちゃめちゃ複雑な曲線でさえも、たった1つの式で表現できる能力があるし、数式は、低い精度から高い精度まで、どんな精度の近似にも当てはまる1つのルールを書き下すことも得意だから、一見すると、どうやって実装したら良いかよくわからない一般式が目に飛び込んできたりする。

しかし、目的が単に曲線を引くことだけであれば、全体の式がわからなくても、簡単な式をたくさんつなぎ合わせればできる。

中学校で習ったとおり、1次関数は直線、2次関数は頂点を境に折り返す放物線、3次関数はさらにもう1つ変曲点を持てる曲線-- という具合に、一般に、式の次数が増えていくと、複雑な線を近似する能力がどんどん上がっていく。

しかし、次数を上げてしまうと、計算するのも大変だし、制御するのも大変だ。

身のまわりでは、絵を描くツール等でよくみかけるベジェ曲線なんかも、すごく綺麗な線だけど、よくつかわれているのは制御しやすい三次関数だ。これをどうつなぎ合わせて作りたい線を引くかということをむしろ問題にしている。

ということからも想像できるように、一度に全体を計算するよりも、簡単な式をつなぎ合わせ、それを伸ばしていくアプローチの方が、プログラムとしては汎用性が高い。 三次ベジェ曲線だけであれば、バーンスタイン多項式というやつ読めなくても実装するのはけっこう簡単。

ベジェ曲線に代表されるように、たくさんの制御点をもとに、曲線を引くアルゴリズムのことを総称してスプライン曲線といい、ベジェ曲線の他にもいくつかある。 スプラインは、普通、簡単な式で表現できるカーブのアルゴリズムと、それらが連続するようにつなげるアルゴリズムの合わせ技でできている。

エルミート曲線 (Hermite Curve)

エルミート曲線は、曲線のアルゴリズムの1つ。 キャットムル-ロム曲線は、このエルミート曲線を連続するようにつなげるアルゴリズムなので、実装するにはまずはこれの理解が必要になる。

「制御点」を指定して多重補完するベジェ曲線とは違い、エルミート曲線は、始点と終点と開始ベクトルと終点ベクトルを入力に与えてあげる。

つまり、始点から始点ベクトルの方向に飛び出していった曲線が、終点付近では、終点ベクトルの向きで向かってきて到着する。そんな関数を考えるのである。

この関数は、始点からひとつめのカーブへ進み、ふたつめのカーブから終点へ進む、という形をしているので、3次多項式(cubic polinominal) で表現できるはずである。

\displaystyle{
f(t) = at^3 + bt^2 + ct + d \ \ \ ( 0 \leq t \leq 1)
}

こんな形をしているはず。 t は、曲線の始まりを0、終わりを1 とした場合の、曲線の進み具合のパラメータ。この式の解が、特定の位置 t における曲線の座標を表している。

式で書くことができたので、あとは4つの係数 a,b,c,d がわかれば、好きなパラメータtを代入してあげるだけで曲線の座標を得ることができる。

というわけで、4つの係数を求めてみる。

始点は t = 0 なので、t に 0を代入すると、その値は始点の座標になる。

\displaystyle{
f(0) = p_0 = d
}

反対に、tに1を代入すると 終点の座標になる。

\displaystyle{
f(1) = p_1 = a + b + c + d
}

位置だけではなく向きのベクトルも与えるので、速度の関数についても考えてみる。 一階微分したこの式が、特定の位置 t における 曲線の傾きの関数になるはずである。

\displaystyle{
f'(t) = 3at^2 + 2bt + c
}

t に 0 を代入すると、始点における曲線の傾きが得られる。

\displaystyle{
f'(0) = v_0 = c
}

t に 1 を代入すると、終点における曲線の傾きになる。

\displaystyle{
f'(1) = v_1 = 3a + 2b + c
}

ここまでの式を使って整理すると、始点&終点の座標と傾きと、係数a,bの関係が以下のようになる。

\displaystyle{
p_1 = a + b + v_0 + p_0
}

\displaystyle{
v_1 = 3a + 2b + v_0
}

あとは、連立方程式を解けば、a,bが求まる。

\displaystyle{
\begin{cases}
a + b = -p_0 + p_1 - v_0 \\
3c_0 + 2c_1 = -v_0 + v_1
\end{cases}
}

\displaystyle{
a = 2p_0 - 2p_1 - v0 + v1  \\
b = -3p_0 + 3p_1 - 2v_0 - v_1
}

これで、t を入力にとるエルミート曲線の関数を、始点/終点それぞれの座標と向きのベクトルから記述することができる。

\displaystyle{
f(t) = (2p_0 - 2p_1 - v0 + v1)t^3 + (-3p_0 + 3p_1 - 2v_0 - v_1)t^2 + v_0t + v_1
}

C#で書くとだいたいこんなかんじ。

    public struct HermitePoly
    {
        // 始点/終点の座標と ベクトルから、曲線の軌道上の座標を返す
        public Vector3 GetPoint(Vector3 p0, Vector3 p1, Vector3 v0, Vector3 v1, float t)
        {
            var c0 = 2f * p0 + -2f * p1 + v0 + v1;
            var c1 = -3f * p0 + 3f * p1 + -2f * v0 - v1;
            var c2 = v0;
            var c3 = p0;

            var t2 = t * t;
            var t3 = t2 * t;
            return c0 * t3 + c1 * t2 + c2 * t + c3;
        }

        // 始点/終点の座標と ベクトルから、曲線の傾きベクトルを返す
        // この関数もつくっておくと3Dの法線とかに使えて便利
        public Vector3 GetTangent(Vector3 p0, Vector3 p1, Vector3 v0, Vector3 v1, float t)
        {
            var c0 = 6f * p0 - 6f * p1 + 3f * v0 + 3f * v1;
            var c1 = -6f * p0 + 6f * p1 - 4f * v0 - 2f * v1;
            var c2 = v0;

            var t2 = t * t;
            return c0 * t2 + c1 * t + c2;
        }
    }

キャットムル-ロム スプライン (Catmull-Rom Spline)

キャットムル-ロム スプラインは、任意の数の制御点を入力して、そこの点をすべて通るような曲線が引けるアルゴリズム。 これは実は、エルミート曲線を連続するようにつなぎ合わせる手法。

エルミート曲線自体は、始点ベクトルと終点ベクトルを自由に入力することができるアルゴリズムだけど、傾きを自由に入力するのではなく、長いスプラインのなかの、前の制御点と次の制御点との傾きを求めて、それを入力に使う。

エルミート曲線が実装できれば、後は簡単。

    public class CatmullRomSpline
    {
        public IReadOnlyList<Vector3> ControlPoints => _controlPoints;

        readonly HermitePoly _poly;
        List<Vector3> _controlPoints = new List<Vector3>
        {
            new Vector3(0f, 0f, 0f),
            new Vector3(1f, 0f, 0f),
            new Vector3(1f, 0f, 1f),
            new Vector3(1f, 0f, 2f),
            // .. 制御点はいくつでも追加できる
        };

        public CatmullRommSpline(HermitePoly poly)
        {
            _poly = poly;
        }

        public Vector3 GetPoint(float t)
        {
            var l = _controlPoints.Count;
            var progress = (l - 1) * t;
            var i = Mathf.FloorToInt(progress);
            var weight = progress - i;

            if (Mathf.Approximately(weight, 0f) && i >= l - 1)
            {
                i = l - 2;
                weight = 1;
            }

            var p0 = _controlPoints[i];
            var p1 = _controlPoints[i + 1];

            Vector3 p2;
            if (i > 0)
            {
                p2 = 0.5f * (ControlPoints[i + 1] - ControlPoints[i - 1]);
            }
            else
            {
                p2 = ControlPoints[i + 1] - ControlPoints[i];
            }

            Vector3 p3;
            if (i < l - 2)
            {
                p3 = 0.5f * (ControlPoints[i + 2] - ControlPoints[i]);
            }
            else
            {
                p3 = ControlPoints[i + 1] - ControlPoints[i];
            }
            return _poly.GetPoint(p0, p1, p2, p3, weight);
        }

これだけで、制御点をすべて通るなめらかな線が引けた。

f:id:hadashia:20171230150523p:plain