数式表示 \LaTeX(MathJax)

最小二乗法による2次ベジェ曲線の近似方法


n個の点xi , yiから二次ベジェ曲線の式を最小二乗法により近似してみます。
具体的には上図の緑の点を始終点とする黒の点付近を通る二次ベジェ曲線である赤線を描画するための青の制御点の座標を近似することです。
各点にばらつきがある場合は誤差がある場合に有効な方法となります。
最小二乗法は読んで字のごとく差の二乗の和が最小となるように係数を計算する方法です。
ただしあくまでも近似ですので、さらに精度を求めるときは最小二乗法の値をもとにさらに別の方法で近似する必要がるかもしれません。
始点をSx , sy、終点をEx , Ey、制御点をCx , Cyとし黒の点の座標をそれぞれxi , yiとします。
x=Ex \times t^2+2 Cx t (1-t)+Sx(1-t)^2
y=Ey \times t^2+2 Cy t (1-t)+Sy(1-t)^2
始終点が既知点となりますので、xi , yiからCx , Cyを最適化したいのですが、媒介変数tが邪魔となります。
最小二乗法により二次方程式のA,B,C,D,Eを求め二次曲線の式から変換可能な場合二次ベジェ曲線に変換します。
x^2+Axy+By^2+Cx+Dy+E=0
上記の式が0にならない場合を誤差とし二乗します。
(By^2+Axy+Dy+x^2+Cx+E)^2
A,B,C,D,Eそれぞれについて偏微分します。
展開すると B^2y^4+2ABxy^3+2BDy^3+2Bx^2y^2+A^2x^2y^2+ 2ADxy^2+2BCxy^2+2BEy^2+D^2y^2+2Ax^3y+2 Dx^2y+2ACx^2y+2AExy+2CDxy+2DEy+x ^4+2Cx^3+2Ex^2+C^2x^2+2CEx+E^2
=B^2y^4+\left(2ABx+2BD\right)y^3+\left(\left(2B+A^2 \right)x^2+\left(2AD+2BC\right)x+2BE+D^2\right)y^2 +\left(2Ax^3+\left(2D+2AC\right)x^2+\left(2AE+2C D\right)x+2DE\right)y+x^4+2Cx^3+\left(2E+C^2\right)x ^2+2CEx+E^2
\displaystyle \frac{\partial }{\partial A}=2Ax^2y^2+2Bxy^3+2Cx^2y+2Dxy^2+2Exy+2x^3y
\displaystyle \frac{\partial }{\partial B}=2Axy^3+2By^4+2Cxy^2+2Dy^3+2Ey^2+2x^2y^2
\displaystyle \frac{\partial }{\partial C}=2Ax^2y+2Bxy^2+2Cx^2+2Dxy+2Ex+2x^3
\displaystyle \frac{\partial }{\partial D}=2Axy^2+2By^3+2Cxy+2Dy^2+2Ey+2x^2y
\displaystyle \frac{\partial }{\partial E}=2Axy+2By^2+2Cx+2Dy+2E+2x^2
行列形式でA,B,C,D,Eについて解く。
A,B,C,D,Eの係数を係数行列Kと置きます。
K \cdot X=U

係数行列K
12345
1Σx2y2 Σxy3 Σx2y Σxy2 Σxy
2Σxy3 Σy4 Σxy2 Σy3 Σy2
3Σx2y Σxy2 Σx2 ΣxyΣx
4Σxy2 Σy3 Σxy Σy2 Σy
5Σxy Σy2 Σx ΣyΣ1

5元連立方程式をここでは逆行列の計算(ピボット選択有)を使用して解きます。
式の定数部を右辺に移動させ定数行列としてUと置きます。

未知数行列X
1
1A
2B
3C
4D
5E

定数行列U
1
1-Σx3y
2-Σx2y2
4-Σx3
3-Σx2y
5-Σx2

X=K^{-1} \cdot U
ax^2+bxy+cy^2+dx+ey+f=0
上記の式で表される線には係数が0になることも考慮すれば、楕円、放物線、双曲線、直線、二重線、点があるので、判別式で判定する。
ax^2+bx+c=0
D=b^2-ac
D>0  異なる2つの実数解
D=0  1つの実数解
D<0  異なる2つの虚数解

直線の判別式

ax^2+(by+d)x+cy^2+ey+f=0
Ax=a
Bx=by+d
Cx=cy^2+ey+f
Dx=(by+d)^2-a(cy^2+ey+f)
cy^2+(bx+e)y+ax^2+dx+f=0
Ay=c
By=(bx+e)
Cy=ax^2+dx+f=0
Dy=(bx+e)^2-c(ax^2+dx+f)
Dx=Dy=0  の時2つの直線を表す。

曲線の判別式

\displaystyle D=B-\frac{A^2}{4}
D > 0楕円
D < 0双曲線
D = 0放物線

放物線

放物線はベジェ曲線に変換可能です。
始終点はそのまま使用し、始終点の接線を求めその交点を制御点とすれば二次ベジェ曲線に置き換えできます。
x^2+Axy+By^2+Cx+Dy+E=0
傾きはxとyの微分値により算出することができます。
dx=Ay+2x+C
dy=2By+Ax+D
x,yを指定して接線の傾きを求めると
\displaystyle d=-\frac{dx}{dy}=-\frac{Ay+2x+C}{2By+Ax+D}
傾きが求まれば切片を算出することができるので直線の式が完成します。
直線の式y=ax+b
始点側の接線
\displaystyle as=-\frac{dx}{dy}=-\frac{Asy+2sx+C}{2Bsy+Asx+D}
bs=sy-as \cdot sx
y=asx+bs
終点側の接線
\displaystyle ae=-\frac{dx}{dy}=-\frac{Aey+2ex+C}{2Bey+Aex+D}
be=ey-es \cdot ex
y=aex+be
交点(cx,cy)
as \cdot x+bs=ae \cdot x+be
(as-ae) x=be-bs
\displaystyle cx=\frac{be-bs}{as-ae}
cy=ae \cdot cx+be

最小二乗法による2次ベジェ曲線の近似の計算例

以下のフォームに座標の一覧を入力して計算ボタンをクリックすると円の近似の結果が表示されます。
入力された各点は黒色の小さい円、近似により算出された円は黒い大きい円、赤色の円が近似された円の中心です。
下表は計算の途中結果を表示しています。一番下の行が合計となります。
座標リスト


係数行列Kの要素の各項の値を算定
i xi yi xiyi xi2 yi2 xi2yi xiyi2 xi2yi2 xi3 yi3 xi3yi xiyi3 xi4 yi4

係数行列K
12345
1Σx2y2 Σxy3 Σx2y Σxy2 Σxy
2Σxy3 Σy4 Σxy2 Σy3 Σy2
3Σx2y Σxy2 Σx2 ΣxyΣx
4Σxy2 Σy3 Σxy Σy2 Σy
5Σxy Σy2 Σx ΣyΣ1
係数行列K
定数行列U
1
1-Σx3y
2-Σx2y2
4-Σx3
3-Σx2y
5-Σx2
定数行列U
係数行列K-1
X=K-1U



始点側の接線




終点側の接線




交点(cx,cy)
as \cdot x+bs=ae \cdot x+be
(as-ae) x=be-bs