3次ベジェ曲線による円弧の近似 - 機械学習基礎理論独習

はじめに

3次ベジェ曲線による円弧の近似はこちら(以下、参考サイトと記載)に必要な事項がすべてまとめられています。
なので、上記のサイトと内容が重複します。

3次ベジェ曲線の制御点を求める

3次ベジェ曲線で円弧を近似します。
円弧は単位円の部分集合とし、円弧の角度を $\theta\in(0,2\pi)$ とします。
両端点 ${\bf P}_0,{\bf P}_3$ は円弧と一致するものとします。
両端点の接線ベクトルは円弧の接線ベクトルと方向が同じものであるとします。
${\bf P}_2-{\bf P}_3$ が $(0,1)$ と同じ方向であることから ${\bf P}_2=(1,\kappa)$ とおくと、 ${\bf P}_1=(\cos\theta+\kappa\sin\theta,\sin\theta-\kappa\cos\theta)$ となることが分かります。

$\begin{eqnarray} {\bf P}(t)=(1-t)^3{\bf P}_0+3(1-t)^2t{\bf P}_1+3(1-t)t^2{\bf P}_2+t^3{\bf P}_3\tag{1} \end{eqnarray}$

$\begin{eqnarray} \left\{ \begin{array}{l} {\bf P}_0=(\cos\theta,\sin\theta)\\ {\bf P}_1=(\cos\theta+\kappa\sin\theta,\sin\theta-\kappa\cos\theta)\\ {\bf P}_2=(1, \kappa)\\ {\bf P}_3=(1,0)\\ \end{array} \right.\tag{2} \end{eqnarray}$

${\bf P}(t)$ の $x,y$ 成分をそれぞれ ${\bf P}_x(t), {\bf P}_y(t)$ とおきます。
曲線上の点 ${\bf P}(1/2)$ は円弧上にあるとします。

$\begin{eqnarray} \left\{\begin{array}{l} {\bf P}_x(1/2)=\cos\dfrac{\theta}{2}\ \ \ \ \ \ (3)\\ {\bf P}_y(1/2)=\sin\dfrac{\theta}{2}\ \ \ \ \ \ (4)\\ \end{array}\right. \end{eqnarray}$

式 $(3)$ を変形していきます。

$\begin{eqnarray} &&{\bf P}_x(1/2)=\cos\frac{\theta}{2}\\ &&\Rightarrow\frac{1}{8}\cos\theta+\frac{3}{8}\left(\cos\theta+\kappa\sin\theta\right)+\frac{3}{8}+\frac{1}{8}=\cos\frac{\theta}{2}\\ &&\Rightarrow4\cos\theta+3\kappa\sin\theta+4=8\cos\frac{\theta}{2}\\ &&\Rightarrow4\left(2\cos^2\frac{\theta}{2}-1\right)+6\kappa\sin\frac{\theta}{2}\cos\frac{\theta}{2}+4=8\cos\frac{\theta}{2}\\ &&\Rightarrow8\cos^2\frac{\theta}{2}+6\kappa\sin\frac{\theta}{2}\cos\frac{\theta}{2}=8\cos\frac{\theta}{2}\\ &&\Rightarrow6\kappa\sin\frac{\theta}{2}\cos\frac{\theta}{2}=8\cos\frac{\theta}{2}-8\cos^2\frac{\theta}{2}\\ &&\Rightarrow\kappa=\frac{4}{3}\left(\frac{\cos\frac{\theta}{2}-\cos^2\frac{\theta}{2}}{\sin\frac{\theta}{2}\cos\frac{\theta}{2}}\right)\ \ \ (\theta\not=\pi)\tag{3.1}\\ &&\Rightarrow\kappa=\frac{4}{3}\left(\frac{1-\cos\frac{\theta}{2}}{\sin\frac{\theta}{2}}\right)\\ &&\Rightarrow\kappa=\frac{4}{3}\left(\frac{1-\left(1-2\sin^2\frac{\theta}{4}\right)}{2\sin\frac{\theta}{4}\cos\frac{\theta}{4}}\right)\\ &&\Rightarrow\kappa=\frac{4}{3}\left(\frac{2\sin^2\frac{\theta}{4}}{2\sin\frac{\theta}{4}\cos\frac{\theta}{4}}\right)\\ &&\Rightarrow\kappa=\frac{4}{3}\left(\frac{\sin\frac{\theta}{4}}{\cos\frac{\theta}{4}}\right)\\ &&\Rightarrow\kappa=\frac{4}{3}\tan\frac{\theta}{4}\tag{3.2}\\ \end{eqnarray}$

式 $(3.1)$ の分母ですが、 $\theta\in(0,2\pi)$ の範囲で $\sin\frac{\theta}{2}\not=0$ ですが、 $\theta=\pi$ で $\cos\frac{\theta}{2} = 0$ となるので、 $\theta\not=\pi$ としてあります。
そこで、 $\theta=\pi$ のときは、式 $(4)$ を変形していきます。

$\begin{eqnarray} &&{\bf P}_y(1/2)=\sin\frac{\pi}{2}\\ &&\Rightarrow\frac{1}{8}\sin\pi+\frac{3}{8}\left(\sin\pi-\kappa\cos\pi\right)+\frac{3}{8}\kappa=1\\ &&\Rightarrow\kappa=\frac{4}{3}\tag{4.1} \end{eqnarray}$

式 $(3.2),(4.1)$ を併せて、以下のように書けます。

$\begin{eqnarray} \kappa=\frac{4}{3}\tan\frac{\theta}{4}\tag{5}\\ \end{eqnarray}$

式 $(5)$ により、 $\kappa$ は $\theta$ にのみ依存することが分かりました。
$\kappa$ が求まったので、3次ベジェ曲線の制御点の座標が定まりました。

円弧の中心からの距離が最大となるパラメータの値

円弧の中心からの距離が最大となるパラメータ $t$ の値を求めてみます。
円弧の中心は原点なので、距離の2乗は以下のようになります。

$\begin{eqnarray} L(t)={\bf P}_x(t)^2+{\bf P}_y(t)^2\tag{6} \end{eqnarray}$

式 $(6)$ を $t$ で微分して $=0$ と置きます。

$\begin{eqnarray} L'(t)=0\tag{7} \end{eqnarray}$

$L'(t)$ はかなり複雑な式なります。(なので、記載しません(できません)。)
式 $(7)$ を $t$ について解いた結果は、以下となります。(この解法については、本記事の下側に記載します。)

$\begin{eqnarray} t=0,\frac{1}{2},1,\frac{1}{2}-\frac{\sqrt{3}}{6},\frac{1}{2}+\frac{\sqrt{3}}{6}\tag{8} \end{eqnarray}$

式 $(8)$ の $t=0,1/2,1$ は円弧上の点であり、円弧の中心からの距離は極小値(最小値) $1$ を取ります。

$\begin{eqnarray} \left|\left| {\bf P}(0) \right|\right|=\left|\left| {\bf P}(1/2) \right|\right|=\left|\left| {\bf P}(1)\right|\right|=1\tag{9} \end{eqnarray}$

$t=\frac{1}{2}\pm\frac{\sqrt{3}}{6}$ の時、円弧の中心からの距離は極大値(最大値)を取ります。

$\begin{eqnarray} \left|\left| {\bf P}\left(\frac{1}{2}\pm\frac{\sqrt{3}}{6}\right)\right|\right|=\sqrt{\frac{10 \sin{\left(\theta \right)} \tan{\left(\frac{\theta}{4} \right)}}{27} - \frac{4 \cos{\left(\theta \right)} \tan^{2}{\left(\frac{\theta}{4} \right)}}{27} + \frac{11 \cos{\left(\theta \right)}}{54} + \frac{8 \tan^{2}{\left(\frac{\theta}{4} \right)}}{27} + \frac{43}{54}}\tag{10} \end{eqnarray}$

式 $(10)$ は参考サイトと違って見えますが、合っていると思います。

$\theta=45,90,135,180,225,270$ 度の時の $\left|\left| {\bf P}(\frac{1}{2}\pm\frac{\sqrt{3}}{6})\right|\right|$ を計算してみます。

角度	長さ
45度	1.00000424552873
90度	1.00027253000743
135度	1.00314575375998
180度	1.01835015443463
225度	1.07638171145720
270度	1.27635594382687

上の表より、3次ベジェ曲線で円弧を近似するときはが90度ぐらいで、区切るのが良さそうです。
グラフは以下になります。

式 (3) のSymPyによる解法

式を立てた後に $\theta \leftarrow 4\alpha$ 式を解いた後に $\alpha\leftarrow\frac{1}{4}\theta$ と戻しているのはそうしないと SymPy では答えがキレイにならないためです。

import sympy

sympy.init_printing()

sympy.var("px p0 p1 p2 p3 t theta kappa alpha")
px = (1-t)**3 * p0 + 3 * (1-t) ** 2 * t * p1 + 3 * (1-t) * t ** 2 * p2 + t ** 3 * p3

px = px.subs([
    (p0, sympy.cos(theta)),
    (p1, sympy.cos(theta) + kappa * sympy.sin(theta)),
    (p2, 1),
    (p3, 1),
    (t, sympy.Rational(1,2))
])

eq = sympy.Eq(px, sympy.cos(theta * sympy.Rational(1,2)))
eq = eq.subs([
    (theta, 4 * alpha)
])
solve = sympy.solve (eq, kappa)
simplified = solve[0].subs([
    (alpha, theta * sympy.Rational(1,4))
])

display(solve)
display(simplified)

コードの実行結果

式 (7) のSymPyによる解法

SymPyでは式 $(7)$ は $\theta$ のままだと $t$ について解けないようです。(参考サイトによるとmaximaなら解けます。)
$\theta$ に適当な値を代入すれば $t$ について解くことができます。

import sympy

sympy.init_printing()

sympy.var("px py p0_x p0_y p1_x p1_y p2_x p2_y p3_x p3_y t theta kappa alpha")
px = (1-t)**3 * p0_x + 3 * (1-t) ** 2 * t * p1_x + 3 * (1-t) * t ** 2 * p2_x + t ** 3 * p3_x
py = (1-t)**3 * p0_y + 3 * (1-t) ** 2 * t * p1_y + 3 * (1-t) * t ** 2 * p2_y + t ** 3 * p3_y

len2 = px ** 2 + py ** 2

len2 = len2.subs([
    (p0_x, sympy.cos(theta)), (p0_y, sympy.sin(theta)),
    (p1_x, sympy.cos(theta) + kappa * sympy.sin(theta)), (p1_y, sympy.sin(theta) - kappa * sympy.cos(theta)),
    (p2_x, 1), (p2_y, kappa),
    (p3_x, 1), (p3_y, 0),
    (kappa, sympy.Rational(4,3) * sympy.tan(theta / 4)),
    (theta, sympy.pi * 0.5) # change here
])

len2d = sympy.diff(len2, t)
len2d = sympy.simplify(len2d)
eq = sympy.Eq(len2d, 0)
solve = sympy.solve(eq, t)

display(solve)

コードの実行結果

式 (10) のSymPyによる解法

import sympy

sympy.init_printing()

sympy.var("px py p0_x p0_y p1_x p1_y p2_x p2_y p3_x p3_y t theta kappa alpha")
px = (1-t)**3 * p0_x + 3 * (1-t) ** 2 * t * p1_x + 3 * (1-t) * t ** 2 * p2_x + t ** 3 * p3_x
py = (1-t)**3 * p0_y + 3 * (1-t) ** 2 * t * p1_y + 3 * (1-t) * t ** 2 * p2_y + t ** 3 * p3_y

len2 = px ** 2 + py ** 2

len2 = len2.subs([
    (p0_x, sympy.cos(theta)), (p0_y, sympy.sin(theta)),
    (p1_x, sympy.cos(theta) + kappa * sympy.sin(theta)), (p1_y, sympy.sin(theta) - kappa * sympy.cos(theta)),
    (p2_x, 1), (p2_y, kappa),
    (p3_x, 1), (p3_y, 0),
    (kappa, sympy.Rational(4,3) * sympy.tan(theta / 4)),
    (t, sympy.Rational(1,2) +(sympy.sqrt(3) / 6))
])

len2 = sympy.simplify(len2)
len = sympy.sqrt(len2)

lensub = len.subs([
    (theta, sympy.pi * 0.5), # change here
])
print(lensub.evalf())
print("\\begin{eqnarray}")
print("&&" + sympy.latex(len2) + "\\\\")
print("\\end{eqnarray}")

コードの実行結果
1.00027253000743
\begin{eqnarray}
&&\frac{10 \sin{\left(\theta \right)} \tan{\left(\frac{\theta}{4} \right)}}{27} - \frac{4 \cos{\left(\theta \right)} \tan^{2}{\left(\frac{\theta}{4} \right)}}{27} + \frac{11 \cos{\left(\theta \right)}}{54} + \frac{8 \tan^{2}{\left(\frac{\theta}{4} \right)}}{27} + \frac{43}{54}\\
\end{eqnarray}

参考リンク

ベジェ曲線で描く円弧と本物の円弧の差（中心角について一般化）