はじめに
Wikipediaの『ベジェ曲線』に、3次ベジェ曲線の変曲点や自己交差の交点の計算式を追記しました。
変曲点にしても交点にしても方程式を解けば求まるので、論文ネタとして魅力がないせいでしょうか、文献がなかなか見当たりません。結局、自力で解くなどしたものを載せたので、導出過程を投稿しておきます。
準備
3次ベジェ曲線を、
\begin{align}
& a_0=-{B_0}_x+3{B_1}_x-3{B_2}_x+{B_3}_x \\
& a_1=3{B_0}_x-6{B_1}_x+3{B_2}_x \\
& a_2=-3{B_0}_x+3{B_1}_x \\
& a_3={B_0}_x \\
& b_0=-{B_0}_y+3{B_1}_y-3{B_2}_y+{B_3}_y \\
& b_1=3{B_0}_y-6{B_1}_y+3{B_2}_y \\
& b_2=-3{B_0}_y+3{B_1}_y \\
& b_3={B_0}_y \\
& x=a_0t^3+a_1t^2+a_2t+a_3 \\
& y=b_0t^3+b_1t^2+b_2t+b_3
\end{align}
とします。また、
\begin{align}
& v=a_0b_1-a_1b_0 \\
& w=a_0b_2-a_2b_0 \\
& u=a_1b_2-a_2b_1
\end{align}
とします。
変曲点
変曲点は、パラメトリック曲線の曲率の定義式
k = \frac{x'y''-y'x''}{\bigl({x'}^2+{y'}^2\bigr)\vphantom{'}^{3/2}}
の分子 $k_a=x'y''-y'x''y$ がゼロになるパラメータ$t_1$と$t_2$です。
$x$と$y$の1次微分と2次微分
\begin{align}
& x'=3a_0t^2+2a_1t+a_2 \\
& x''=6a_0t+2a_1 \\
& y'=3b_0t^2+2b_1t+b_2 \\
& y''=6b_0t+2b_1
\end{align}
を$ka$に代入すると、
\begin{align}
ka & = -2(3b_0t^2+2b_1t+b_2)(3a_0t+a_1)+2(3a_0t^2+2a_1t+a_2)(3b_0t+b_1) \\
& = (6a_1b_0-6a_0b_1)t^2+(6a_2b_0-6a_0b_2)t+2a_2b_1-2a_1b_2 \\
& = 3(a_1b_0-a_0b_1)t^2+3(a_2b_0-a_0b_2)t+a_2b_1-a_1b_2 \\
& = 0
\end{align}
$v, w, u$ を代入します。
\begin{align}
& -3vt^2-3wt-u=0 \\
& 3vt^2+3wt+u=0
\end{align}
2次方程式の解の公式にあてはめます。
\begin{align}
& t_1=\frac{-b + \sqrt{b^2-4ac}}{2a} \\
& t_2=\frac{-b - \sqrt{b^2-4ac}}{2a} \\
& a=3v \\
& b=3w \\
& c=u \\
& t_1=\frac{-3w + \sqrt{9w^2-12uv}}{6v} \\
& t_2=\frac{-3w - \sqrt{9w^2-12uv}}{6v}
\end{align}
$t_1\in[0,1]$を満たす場合、$t_1$は変曲点です。
$t_2\in[0,1]$を満たす場合、$t_2$は変曲点です。
（変曲点は0～2個です）
var ('a0 a1 a2 a3 b0 b1 b2 b3')
v=a0*b1-a1*b0
w=a0*b2-a2*b0
u=a1*b2-a2*b1
t1=(-3*w + sqrt(9*w^2 - 12*u*v))/(6*v)
t2=(-3*w - sqrt(9*w^2 - 12*u*v))/(6*v)
t=t1
x=a0*t^3+a1*t^2+a2*t+a3
y=b0*t^3+b1*t^2+b2*t+b3
x1=3*a0*t^2+2*a1*t+a2
y1=3*b0*t^2+2*b1*t+b2
x2=6*a0*t+2*a1
y2=6*b0*t+2*b1
ka1=(x1*y2-y1*x2)
t=t2
x=a0*t^3+a1*t^2+a2*t+a3
y=b0*t^3+b1*t^2+b2*t+b3
x1=3*a0*t^2+2*a1*t+a2
y1=3*b0*t^2+2*b1*t+b2
x2=6*a0*t+2*a1
y2=6*b0*t+2*b1
ka2=(x1*y2-y1*x2)
ka1.simplify_full()
ka2.simplify_full()
0
0
自己交差の交点
交点は、自力で解こうとするとなかなか厄介で、範囲を広げて探したところ、以下の2件がありました。
-
Stack Exchange - Mathematics
2D cubic Bezier curve. Point of self-intersection
https://math.stackexchange.com/questions/3776840/2d-cubic-bezier-curve-point-of-self-intersection -
Zenn 白山風露
自己交差する3次ベジェ曲線の交差点を求める
https://zenn.dev/kazatsuyu/scraps/a362d9d1c35da6
1)は最後まで書かれておらず、2)を全面的に参考にさせていただきました。
ただし、誤差が生じやすい演算（高次項どうしの減算）が含まれていたため、修正したものをWikipediaに載せました。
まず、白山さんの解法そのままです。ただし変数名などはWikipediaの記事に合わせて変更してあります。
3次ベジェ曲線 $x=f(t), y=g(t)$ にループがある場合の交点は、
$f(t_1)=f(t_2), g(t_1)=g(t_2)$となるパラメータ$t_1$と$t_2$とします。
$x$については、
$f(t)=a_0 t^3+a_1 t^2+a_2 t+a_3, f(t_1)=f(t_2)$
なので
$a_0 t_1^3+a_1 t_1^2+a_2 t_1+a_3=a_0 t_2^3+a_1 t_2^2+a_2 t_2+a_3$
左辺に集めると
$a_0 t_1^3+a_1 t_1^2+a_2 t_1+a_3-a_0 t_2^3-a_1 t_2^2-a_2 t_2-a_3=0$
$a_0 (t_1^3-t_2^3)+a_1 (t_1^2-t_2^2)+a_2 (t_1-t_2)=0$
ここで
$(t_1^3-t_2^3)=(t_1^2+t_1 t_2+t_2^2) (t_1-t_2)$
$(t_1^2-t_2^2)=(t_1+t_2) (t_1-t_2)$
なので、両辺を$(t_1-t_2)$で割ると
$a_0 (t_1^2+t_1 t_2+t_2^2)+a_1 (t_1+t_2)+a_2=0$ (1)
$y$についても同様に
$b_0 (t_1^2+t_1 t_2+t_2^2)+b_1 (t_1+t_2)+b_2=0$ (2)
(1)を変形します。
$t_1^2+t_1 t_2+t_2^2 = -(a_1 (t_1+t_2)+a_2)/a_0$
(2)に代入します。
$b_0/a_0 (-a_1 (t_1+t_2)-a_2)+b_1 (t_1+t_2)+b_2=0$
$t_1$について整理します。
$-b_0 a_1 t_1/a_0-b_0 a_1 t_2/a_0-b_0 a_2/a_0+b_1 t_1+b_1 t_2+b_2=0$
$-(b_0 a_1/a_0-b_1) t_1=b_0 a_1 t_2/a_0+b_0 a_2/a_0-b_1 t_2-b_2$
$-(b_0 a_1/a_0-b_1) t_1=(b_0 a_1/a_0-b_1) t_2+b_0 a_2/a_0-b_2$
$t_1=-t_2-(b_0 a_2/a_0-b_2)/(b_0 a_1/a_0-b_1)$
$t_1=-t_2-(b_0 a_2-b_2 a_0)/(b_0 a_1-b_1 a_0)$
$t_1=-t_2-(-w)/(-v)$
$t_1=-t_2-w/v$ (3)
(1)に代入します。
$a_0 ((-t_2-w/v)^2+(-t_2-w/v) t_2+t_2^2)+a_1 ((-t_2-w/v)+t_2)+a_2=0$
$a_0 (t_2^2+2 t_2 w/v+w^2/v^2-t_2^2-w t_2/v+t_2^2)+a_1 (-t_2-w/v+t_2)+a_2=0$
$a_0 t_2^2-a_0 t_2^2+a_0 t_2^2 +a_0 2 t_2 w/v-a_0 w t_2/v +a_0 w^2/v^2 -a_1 w/v+a_2=0$
$(a_0-a_0+a_0) t_2^2 +(2 a_0 w/v-a_0 w/v) t_2 +a_0 w^2/v^2 -a_1 w/v+a_2=0$
$a_0 t_2^2 + a_0 w t_2/v + a_0 w^2/v^2-a_1 w/v+a_2=0$
両辺に$v^2$を掛けます。
$a_0 v^2 t_2^2 + a_0 w v^2 t_2/v + a_0 w^2 v^2/v^2-a_1 w v^2/v+a_2 v^2=0$
$a_0 v^2 t_2^2 + a_0 w v t_2 + a_0 w^2 -a_1 w v +a_2 v^2=0$
2次方程式の解の公式にあてはめます。
\begin{align}
& t_2 =\frac{-b \pm \sqrt{b^2-4ac}}{2a} \\
& a =a_0 v^2 \\
& b =a_0 v w \\
& c =a_0 w^2-a_1 v w+a_2 v^2 \\
\end{align}
(3)に代入します。
\begin{align}
t_1 & =-\frac{-b \pm \sqrt{b^2 - 4 a c}}{2a}-\frac{w}{v} \\
& =\frac{b}{2a} \mp \frac{\sqrt{b^2 - 4 a c}}{2a} - \frac{w}{v} \\
& =\frac{a_0 v w}{2 a_0 v^2} - \frac{w}{v} \mp \frac{\sqrt{b^2 - 4 a c}}{2a} \\
& =\frac{w}{2 v} - \frac{w}{v} \mp \frac{\sqrt{b^2 - 4 a c}}{2a} \\
& =-\frac{w}{2 v} \mp \frac{\sqrt{b^2 - 4 a c}}{2a}
\end{align}
$b/a=w/v$なので
\begin{align}
t_1 & =-\frac{b}{2a} \mp \frac{\sqrt{b^2 - 4 a c}}{2a} \\
& =\frac{-b \mp \sqrt{b^2 - 4 a c}}{2a}
\end{align}
var ('a0 a1 a2 a3 b0 b1 b2 b3')
v=a0*b1-a1*b0
w=a0*b2-a2*b0
a=a0*v^2
b=a0*v*w
c=a0*w^2-a1*v*w+a2*v^2
t1=(-b+sqrt(b^2-4*a*c))/2/a
t2=(-b-sqrt(b^2-4*a*c))/2/a
x1=a0*t1^3+a1*t1^2+a2*t1+a3
y1=b0*t1^3+b1*t1^2+b2*t1+b3
x2=a0*t2^3+a1*t2^2+a2*t2+a3
y2=b0*t2^3+b1*t2^2+b2*t2+b3
x=x1-x2
y=y1-y2
x.simplify_full()
y.simplify_full()
0
0
ところが、実際に制御点の座標値をあてはめた検証で問題が発生しました。
例として、以下はWikipediaに載せた類型の内、カスプの制御点の座標値です。
カスプの場合は$t_1=t_2$となることが期待されます。
| 制御点 | $x$ | $y$ |
|---|---|---|
| $B_0$ | $15.45$ | $10.44$ |
| $B_1$ | $8.77$ | $2.13$ |
| $B_2$ | $16.91$ | $2.18$ |
| $B_3$ | $7.155375845$ | $11.44$ |
これをあてはめて交点を計算すると、$t_1=0.485048947, t_2=0.485048947$ となり、問題ありません。
ところが、$x$と$y$を入れ替えたところ、$t_1=0.485049043, t_2=0.485048850$ となりました。$t_1=t_2$とは言い難い差です。
そこで各変数の値を見ると、$b^2-4ac$で計算誤差が発生していました。
| 項 | 値 |
|---|---|
| $b^2$ | $368954101702.591$ |
| $4ac$ | $368954101702.576$ |
| $b^2-4ac$ | $0.014587402$ |
$b^2$と$4ac$を式展開すると
\begin{align}
b^2 & =(a_0 v w)^2 \\
& =(a_0 (a_0 b_1-a_1 b_0) (a_0 b_2-a_2 b_0))^2 \\
4 a c & =4 (a_0 v^2) (a_0 w^2-a_1 v w+a_2 v^2) \\
& =4 (a_0 (a_0 b_1-a_1 b_0)^2) (a_0 (a_0 b_2-a_2 b_0)^2-a_1 (a_0 b_1-a_1 b_0) (a_0 b_2-a_2 b_0)+a_2 (a_0 b_1-a_1 b_0)^2)
\end{align}
となり、座標値の10乗どうしの引算で、64ビット倍精度なのに有効桁が2桁しか残りませんでした。
計算誤差低減版
(3)を(1)に代入するところまで同じです。
$a_0 t_2^2 + a_0 w t_2/v + a_0 w^2/v^2-a_1 w/v+a_2=0$
両辺に（$v^2$ではなく）$v/a_0$を掛けます。
$v t_2^2 + w t_2 + w^2 /v - a_1 w/a_0 + a_2 v/a_0=0$
2次方程式の解の公式にあてはめます。
\begin{align}
t_2 & =\frac{-b \pm \sqrt{b^2-4ac}}{2a} \\
a & =v \\
b & =w \\
c & =w^2/v-a_1 w/a_0+a_2 v/a_0 \\
t_2 & =\frac{-w \pm \sqrt{w^2 - 4 v (w^2/v - a_1 w/a_0 + a_2 v/a_0)}}{2 v} \\
& =\frac{-w \pm \sqrt{w^2 - 4 w^2 + 4 v a_1 w/a_0 - 4 v a_2 v/a_0}}{2 v} \\
& =\frac{-w \pm \sqrt{-3 w^2 - 4 v (-a_1 w/a_0 + a_2 v/a_0)}}{2 v} \\
& =\frac{-w \pm \sqrt{-3 w^2 - 4 v (-a_1 (a_0 b_2-a_2 b_0)/a_0 + a_2 (a_0 b_1-a_1 b_0)/a_0)}}{2 v} \\
& =\frac{-w \pm \sqrt{-3 w^2 - 4 v (-a_1 a_0 b_2/a_0 +a_1 a_2 b_0/a_0 + a_0 b_1 a_2/a_0 -a_1 b_0 a_2/a_0)}}{2 v} \\
& =\frac{-w \pm \sqrt{-3 w^2 - 4 v (-a_1 b_2 +a_1 a_2 b_0/a_0 + b_1 a_2 -a_1 b_0 a_2/a_0)}}{2 v} \\
& =\frac{-w \pm \sqrt{-3 w^2 + 4 v u}}{2 v}
\end{align}
(3)に代入します。
\begin{align}
& t_1 =-\frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}- \frac{w}{v} \\
\end{align}
$a=v, b=w$なので
\begin{align}
t_1 & =-\frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}- \frac{b}{a} \\
& =-\frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}- \frac{2 b}{2 a} \\
& =\frac{b - 2 b \mp \sqrt{b^2 - 4 a c}}{2 a} \\
& =\frac{-b \mp \sqrt{b^2 - 4 a c}}{2 a}
\end{align}
var ('a0 a1 a2 a3 b0 b1 b2 b3')
v=a0*b1-a1*b0
w=a0*b2-a2*b0
u=a1*b2-a2*b1
t1=1/2*(-w + sqrt(4*u*v - 3*w^2))/v
t2=1/2*(-w - sqrt(4*u*v - 3*w^2))/v
t=t1
x1=a0*t^3+a1*t^2+a2*t+a3
y1=b0*t^3+b1*t^2+b2*t+b3
t=t2
x2=a0*t^3+a1*t^2+a2*t+a3
y2=b0*t^3+b1*t^2+b2*t+b3
x=x1-x2
y=y1-y2
x.simplify_full()
y.simplify_full()
0
0
先ほどと同じ座標値をあてはめて交点を計算すると、$t_1=0.485048947, t_2=0.485048947$ となり、$x$と$y$を入れ替えても同じく、$t_1=0.485048947, t_2=0.485048947$ となりました。
| 項 | 値 |
|---|---|
| $4uv$ | $2079716.139$ |
| $3w^2$ | $2079716.139$ |
| $4uv-3w^2$ | $0.000000000$ |
$4uv$と$3w^2$を式展開すると
\begin{align}
& 4uv = 4(a_1b_2-a_2b_1)(a_0b_1-a_1b_0) \\
& 3w^2 =3(a_0b_2-a_2b_0)^2
\end{align}
で、座標値の4乗どうしの引算となり、誤差が低減されました。