【C++】2次のベジエ曲線の補間

C++

はじめに

先日、sablib向けにC++でベジエ曲線のコードを書く機会がありました。そこでは、いくつかの制御点が与えられた状況下で「$x$が等間隔になるよう補間する」ことが必要でしたので、そのあたりの実装方法について少しまとめてみました。

ベジエ曲線とは

ベジエ曲線は、いくつかの「制御点(control points)」によって形が決まる曲線です。曲線自体は始点と終点を必ず通りますが、途中の制御点は曲線の形を引っ張るように影響を与えるだけで、曲線がその点を通過するわけではありません。

制御点の数によって、以下のように分類されます。

  1. 1次ベジエ曲線
    制御点がなく始点と終点のみで、単純な直線になります。
  2. 2次ベジエ曲線
    始点と終点のほかに制御点が1つで、なめらかな1つのカーブを描きます。放物線のような形になります。
  3. 3次ベジエ曲線
    始点と終点のほかに制御点が2つで、S字型を含む、より複雑な曲線を表現できます。

グラフィック関連ではよく3次のベジエ曲線が使われますが、ここではスペクトル等のバックグラウンド推定で使うことを想定して2次のベジエ曲線を考えることとします。

数式による表現

$n$次のベジエ曲線は、以下のバーンスタイン多項式(Bernstein polynomial)を使って表されます。

$$
\displaystyle
B(t) = \sum_{i=0}^{n} \binom{n}{i} (1-t)^{n-i} t^{i} P_{i} \qquad (0 \leqq t \leqq 1)
$$

ここで$P_{i}$は制御点、$t$は0から1まで変化するパラメータです。$t$を少しずつ増やしながら計算した点をつなげていくことで、曲線が描かれます。

2次ベジエ曲線の場合、具体的には次のようになります。

$$
\displaystyle
B(t) = (1-t)^{2} P_{0} + 2(1-t) t P_{1} + t^{2} P_{2}
$$

複数の制御点をつなげて曲線を描く方法

開始点$(x_{s}, y_{s})$と終了点$(x_{e}, y_{e})$、その間にある$m$個の制御点の組$(x_{0}, y_{0}), (x_{1}, y_{1}), \dots, (x_{m}, y_{m})$が与えられたとき、これを1つのベジエ曲線で描こうとすると、高次のベジエ曲線となり計算が複雑になります。そこで実用上よく使われるのが、「曲線をいくつかの2次ベジエ曲線に分割してつなげる」という方法です。

具体的には、以下のような方法を取ります。

1. 制御点の間に終点兼始点となる中間点を挿入する

隣り合う制御点$(x_{i}, y_{i})$と$(x_{i+1}, y_{i+1})$の中点を取ります。

$$
\displaystyle
M_{i} = \left(\frac{x_{i} + x_{i+1}}{2}, \frac{y_{i} + y_{i+1}}{2}\right)
$$

これを、すべての制御点の隣接ペアの間に挿入します。この中間点$M_{i}$は、直前のセグメントの「終点」であると同時に、次のセグメントの「始点」を兼ねる役割を持ちます。制御点同士を直接つなぐのではなく、この中間点を経由させることで、セグメントの境目で曲線の傾きが連続的になり、折れ目のないなめらかな曲線を実現できます。

セグメントに分割した2次のベジエ曲線のイメージ

2. 挿入した中間点を利用して、始点−制御点−終点の組を作り、2次のベジエ曲線を描く

中間点を挿入し終えたら、次のように「始点・制御点・終点」の3点セットをセグメントごとに作っていきます。

  • 1番目のセグメント:開始点$(x_{s}, y_{s})$、制御点$(x_{0}, y_{0})$、中間点$M_{0}$
  • 2番目のセグメント:開始点$M_{0}$、制御点$(x_{1}, y_{1})$、中間点$M_{1}$
  • $\dots$(制御点の数だけ繰り返す)$\dots$
  • $m$番目のセグメント:開始点$M_{m-1}$、制御点$(x_{m}, y_{m})$、終了点$(x_{e}, y_{e})$

こうして作られた3点セットそれぞれに対して、先ほどの2次ベジエ曲線の式を適用して描画すれば、$m$個の制御点をすべて反映しながら、全体として滑らかにつながった1本の曲線が完成します。

$x$座標を等間隔にして補間するには

バーンスタイン多項式は$x$座標と$y$座標に分割することもできます。

$$
\displaystyle
x(t) = (1-t)^{2} x_{0} + 2(1-t) t x_{1} + t^{2} x_{2} \\
y(t) = (1-t)^{2} y_{0} + 2(1-t) t y_{1} + t^{2} y_{2}
$$

ここで注意したいのは、$t$を$0, 0.1, 0.2, \dots, 1.0$のように均等に増やしても、それに対応する$x(t)$の値は均等な間隔にはならないという点です。$x(t)$は$t$の2次関数であり、直線的な対応関係にないためです。

このような場合、「欲しい$x$の値に対応する$t$」を先に求めてから、その$t$を使って$y(t)$を計算する、という順番の処理が必要になります。

スペクトルやクロマトグラムのようなデータであれば、$x$について等間隔に並んでいるため、補間をする際に少し困ったことになります。

2次ベジエ曲線であれば、$x(t)$の式を$t$について整理すると

$$
x(t) = a t^2 + b t + c
$$

$$
a = x_{0} – 2x_{1} + x_{2}, \quad b = 2(x_{1} – x_{0}), \quad c = x_{0}
$$

という単純な2次方程式の形になります。求めたい$x$の値を左辺に代入すれば、あとは2次方程式の解の公式

$$
t = \frac{-b + \sqrt{b^2 – 4a(c – x)}}{2a}
$$

を使って直接$t$を求めることができます($a=0$になる場合は1次方程式として別途処理します)。この$t$を$y(t)$の式に代入すれば、目的の$x$に対応する$y$が得られます。

3次ベジエ曲線になると $x(t)$は3次方程式になり、解の公式は一気に複雑になります。そこで実用上は、ニュートン法などを使って数値的に$t$を求める方法がよく使われます。

C++による実装例

いくつかの制御点が与えられた状況下で、2次のベジエ曲線を作成し$x$が等間隔となるよう補間するコードの実装例を以下に示します。

#include <array>
#include <cmath>
#include <iostream>
#include <vector>
#include <algorithm>

struct Point {
    double x;
    double y;
};

// 2次ベジエ曲線
Point Bezier(const Point& p0, const Point& p1, const Point& p2, double t)
{
    const double u = 1.0 - t;

    return {
        u * u * p0.x + 2.0 * u * t * p1.x + t * t * p2.x,
        u * u * p0.y + 2.0 * u * t * p1.y + t * t * p2.y
    };
}

// xからtを求める
double Solve(const Point& p0, const Point& p1, const Point& p2, double x)
{
    const double a = p0.x - 2.0 * p1.x + p2.x;
    const double b = 2.0 * (p1.x - p0.x);
    const double c = p0.x - x;

    // 直線の場合
    if (std::abs(a) < 1e-12) {
        return -c / b;
    }

    const double d = b * b - 4.0 * a * c;

    if (d < 0.0) {
        return 0.0;
    }

    const double s = std::sqrt(d);

    const double t1 = (-b + s) / (2.0 * a);
    const double t2 = (-b - s) / (2.0 * a);

    // 通常は [0, 1] に入る解を使用
    if (0.0 <= t1 && t1 <= 1.0) {
        return t1;
    }

    return t2;
}

// 制御点を使って滑らかな2次ベジエ曲線を構築し、
// x が等間隔になるよう補間する
std::vector<Point> BezierInterpolate(
    const Point& start, const std::vector<Point>& controls,
    const Point& end, std::size_t n
)
{
    std::vector<Point> result;

    if (n == 0) {
        return result;
    }

    // 中間点を作る
    std::vector<Point> mid;

    for (std::size_t i = 0; i + 1 < controls.size(); i++) {
        mid.push_back({
            (controls[i].x + controls[i + 1].x) / 2.0,
            (controls[i].y + controls[i + 1].y) / 2.0
        });
    }

    // 各2次ベジエセグメント
    std::vector<std::array<Point, 3>> segments;

    for (std::size_t i = 0; i < controls.size(); i++) {
        Point p0;
        Point p2;

        if (i == 0) {
            p0 = start;
        }
        else {
            p0 = mid[i - 1];
        }

        if (i + 1 == controls.size()) {
            p2 = end;
        }
        else {
            p2 = mid[i];
        }

        segments.push_back({p0, controls[i], p2});
    }

    // x の最小値~最大値を等間隔にサンプリング
    const double xmin = start.x;
    const double xmax = end.x;

    for (std::size_t i = 0; i < n; i++) {
        const double x = xmin + (xmax - xmin) * i / (n - 1);

        // x を含むセグメントを探す
        for (const auto& seg : segments) {
            if (x >= seg[0].x && x <= seg[2].x) {
                const double t = Solve(seg[0], seg[1], seg[2], x);
                result.push_back(Bezier(seg[0], seg[1], seg[2], t));
                break;
            }
        }
    }

    return result;
}

int main()
{
    Point start{0.0, 1.0};

    std::vector<Point> controls{
        {1.0, 3.0},
        {2.0, 0.5},
        {4.0, 2.5}
    };

    Point end{5.0, 1.5};

    // x方向に100点補間
    const auto result = BezierInterpolate(start, controls, end, 100);

    for (const auto& p : result) {
        std::cout << p.x << ", " << p.y << '\n';
    }
}

まとめ

今回は、複数の制御点から滑らかな曲線を作成する方法として、2次ベジエ曲線を使った補間について説明しました。

複数の制御点をそのまま1本のベジエ曲線として扱うのではなく、制御点の間に中間点を挿入して複数の2次ベジエ曲線に分割することで、比較的簡単に滑らかな曲線を構成できます。

また、スペクトルやクロマトグラムのように$x$が等間隔に並んだデータへ適用する場合、単純にパラメータ$t$を等間隔に変化させるだけでは$x$も等間隔にはなりません。そのため、まず目的の$x$に対応する$t$を求め、その$t$から$y$を計算する必要があります。

2次ベジエ曲線では$x(t)$が2次方程式になるため、解の公式を利用して$t$を直接求めることができます。一方、3次以上のベジエ曲線ではより複雑になるため、ニュートン法などの数値計算法を利用するのがよいです。

このように、ベジエ曲線は単に曲線を描くだけでなく、スペクトルやクロマトグラムのようなデータの補間にも利用できます。とくに2次ベジエ曲線であれば計算も比較的単純なので、バックグラウンド推定などへの応用にも扱いやすい方法だと思います。

C++

Posted by izadori