ベジェ曲線で円弧を描く
コンピュータグラフィックスで滑らかな曲線を描くのに利用されるベジェ曲線（Bézier Curve）を用いて円（1/4円弧）を描きます。3次ベジェ曲線を用いるので制御点は4つになります。「ベジェ曲線の定義と4つの性質（高校数学の美しい物語）」 で紹介されている以下の２つを用いると4つの制御点は$[1,0],[1,v],[v,1],[0,1]$と表せることが分かるのでこの$v$を求めることを考えます。
- ベジェ曲線の端点
- ベジェ曲線の端点で進む向き
ベジェ曲線の式
ベジェ曲線は一般には以下のように表され、今回は$N=3$のケースです。
\vec{p}(t) = \sum_{k=0}^{N}\binom{N}{k}t^k(1-t)^{N-k}\vec{p_k}
これをそのままpythonのコードにします。bzr_xyは制御点のリストを渡すと$[x,y]$を$t$の式として返します。下の例では$v=0.5$として4点を渡して$[x,y]$の式を$t$を変数として得ています。この式を元にして$t=[0..1]$で描いたのが上の図の緑の曲線です。青の線で描かれた円弧より少し内側になっているのが分かります。
import sympy as sp
import numpy as np
t, v = sp.symbols('t v')
def bzr_xy(p):
N = len(p)-1
return sum([sp.binomial(N,i)*(1-t)**(N-i)*t**i*p[i] for i in range(N+1)])
p = np.array([[1,0],[1,0.5],[0.5,1],[0,1]])
print(bzr_xy(p))
# [0.5*t**2*(3 - 3*t) + 3.0*t*(1 - t)**2 + 1.0*(1 - t)**3
# 1.0*t**3 + 1.0*t**2*(3 - 3*t) + 1.5*t*(1 - t)**2]
45度の半径が1となるようにvを求める
45度の直線$y=x$と円弧の交点は$(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})$となるはずなのでこれを条件として$v$を求めます。$v = 0.552284749830793$が求まりました。
x, y = bzr_xy(np.array([[1,0],[1,v],[v,1],[0,1]]))
sol = sp.solve((x - sp.sqrt(2)/2, y - sp.sqrt(2)/2), (t,v))
v_value = sol[0][1].evalf()
print("v = ", v_value)
# v = 0.552284749830793
ベジェ曲線の円の面積
この時の円の面積を求めます。$1/4$円の面積$A$は以下のように表されます。
A = \int_{0}^{1} y \,dx = \int_{1}^{0} y \frac{dx}{dt} \,dt
これをそのままsympyで書くと以下のようになります。面積の誤差は0.028%となりほぼ等しいことが分かります。
A = sp.integrate(y*sp.diff(x,t), (t, 1, 0))
A_value = A.subs(v,v_value)
print("A = ", A_value, ((A_value*4-sp.pi)*100/sp.pi).evalf())
# A = 0.785618083164127 0.0280010543603986
円の面積と等しくなる時のvを求める
さらに円の面積と等しくなる時の$v$をsympyを使って求めてみます。$A-\pi/4$=0として$v$を求めます。$v = 0.551778477804468$が求まりました。
eq = A - sp.pi/4
vs = sp.solve(eq,v)
print(vs[0].evalf())
# 0.551778477804468
（開発環境：Google Colab）
この考え方はProject Euler、Problem 363: Bézier Curvesを解くのに役に立ちます。