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

n個の点x_i , y_iから1本の二次ベジェ曲線の式を最小二乗法により近似してみます。(1本と断っているのは凹凸が複数ある曲線は1本の二次ベジェ曲線で表現できないからです。)
具体的には上図の緑の点を始終点とする黒の点付近を通る二次ベジェ曲線である赤線を描画するための青の制御点の座標を近似することです。
各点にばらつきがある場合は誤差がある場合に有効な方法となります。
最小二乗法は読んで字のごとく差の二乗の和が最小となるように係数を計算する方法です。
ただしあくまでも近似ですので、さらに精度を求めるときは最小二乗法の値をもとにさらに別の方法で近似する必要がるかもしれません。
始点をSx , sy、終点をEx , Ey、制御点をCx , Cyとし黒の点の座標をそれぞれx_i , y_iとします。
$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
	1	2	3	4	5
1	Σx²y²	Σxy³	Σx²y	Σxy²	Σxy
2	Σxy³	Σy⁴	Σxy²	Σy³	Σy²
3	Σx²y	Σxy²	Σx²	Σxy	Σx
4	Σxy²	Σy³	Σxy	Σy²	Σy
5	Σxy	Σy²	Σx	Σy	Σ1

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

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

定数行列U
	1
1	-Σx³y
2	-Σx²y²
4	-Σx³
3	-Σx²y
5	-Σx²

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

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

係数行列Kの要素の各項の値を算定
i	x_i	y_i	x_iy_i	x_i²	y_i²	x_i²y_i	x_iy_i²	x_i²y_i²	x_i³	y_i³	x_i³y_i	x_iy_i³	x_i⁴	y_i⁴

係数行列K
	1	2	3	4	5
1	Σx²y²	Σxy³	Σx²y	Σxy²	Σxy
2	Σxy³	Σy⁴	Σxy²	Σy³	Σy²
3	Σx²y	Σxy²	Σx²	Σxy	Σx
4	Σxy²	Σy³	Σxy	Σy²	Σy
5	Σxy	Σy²	Σx	Σy	Σ1

定数行列U
	1
1	-Σx³y
2	-Σx²y²
4	-Σx³
3	-Σx²y
5	-Σx²

始点側の接線

終点側の接線

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

グラフの赤丸が始終点、緑丸が既知点、赤線が近似されたベジェ曲線、青線が制御点と始終点を結んだ線となります。

入力フォームへ移動

以下に近似されたベジェ曲線と既知点との誤差を表に示しています。
t列が既知点に一番近い近似されたベジェ曲線の点(bx,by)の媒介変数を示します。
dがベジェ曲線の点(bx,by)と既知点(xi,yi)との距離すなわち誤差を示します。
表の後には最大の誤差の値を算出し示しています。この値が実用上許容できる値の範囲に収まっていればベジェ曲線の近似に成功していることを示します。

既知点からの距離
i	x_i	y_i	t	bx	by	d

最大の誤差
ベジェ曲線の長さ(SVGのpathの長さをgetTotalLength()で取得した値です。ブラウザにより多少異なった値が取得されます。)
折れ線の長さ

うまくいかない近似例

下表のxi,yiの様に既知点が一直線上に並んでいる場合、上記の方法で近似すると下表の下の図のように制御点が極端な位置で近似されてしまうケースがあります。
ベジェ曲線は媒介変数tが0～1の間で均等に介在するべきですが、ほぼ下表のtで示す通り0～0.1の間にエネルギーが集中しています。

既知点からの距離
i	x_i	y_i	t	bx	by	d
1	4237.0000	5327.0000	2.2204460e-16	4237.0000	5327.0000	1.2862197e-12
2	4260.0000	5337.0000	0.0049509460	4254.3045	5343.1157	8.3570896
3	4281.0000	5350.0000	0.010070127	4272.0346	5359.6272	13.155281
4	4302.0000	5365.0000	0.015529972	4290.7627	5377.0674	16.489332
5	4321.0000	5382.0000	0.021023723	4309.4176	5394.4386	16.996197
6	4341.0000	5402.0000	0.027182531	4330.1047	5413.7013	15.988366
7	4361.0000	5424.0000	0.033716199	4351.7897	5433.8921	13.516053
8	4382.0000	5448.0000	0.040803560	4375.0084	5455.5097	10.260529
9	4404.0000	5473.0000	0.048312213	4399.2620	5478.0894	6.9534855
10	4428.0000	5500.0000	0.056588981	4425.5852	5502.5941	3.5441459
11	4455.0000	5528.0000	0.065684698	4454.0152	5529.0580	1.4454296
12	4485.0000	5557.0000	0.075645895	4484.5525	5557.4808	0.65687467
13	4518.0000	5588.0000	0.086695866	4517.6962	5588.3264	0.44588397
14	4556.0000	5619.0000	0.098939356	4553.5216	5621.6636	3.6382684
15	4598.0000	5651.0000	1.0000000	4598.0000	5651.0000	0.0000000

近似された曲線始点 4237,5327 制御点 5992.4,6961.8 終点 4598,5651
ベジェ曲線の長さ2183.58154296875(SVGのpathの長さをgetTotalLength()で取得した値です。ブラウザにより多少異なった値が取得されます。)
折れ線の長さ488.15524
下図の赤丸が始終点、緑丸がx_i,y_i赤線が近似されたベジェ曲線です。

うまくいかない近似例

上記の様なケースを検出するためにfx(t) fy(t)の極限値をもとめその値が終点の1個手前のtと比較して極限値の時のtの方が大きい場合、うまくいっていないと判断します。
極限値はfx(t)及びfy(t)をそれぞれ微分しfx'(t)=0 fy'(t)=0の時のtを解けば求まります。
$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$
$fx'(t)=2Ex \times t-2Cx \times t-2Sx(1-t)+2Cx(1-t)$
$fy'(t)=2Ey \times t-2Cy \times t-2Sy(1-t)+2Cy(1-t)$
$\displaystyle tx=\frac{Sx-Cx}{Sx+Ex-2Cx}$
$\displaystyle ty=\frac{Sy-Cy}{Sy+Ey-2Cy}$
上式のtx,tyがxとyのそれぞれの極限値となります。
極限値が0～1に収まっていない場合はベジェ曲線の近似に失敗していないことにしました。
誤差が最大の点の倍の位置でベジェ曲線の近似範囲を分割します。
2分割して2つの2次曲線で近似した場合下図の様になります。

係数行列Kの要素の各項の値を算定
i	x_i	y_i	誤差	図形
1	4237.0000	5327.0000	0.0000000	Q4306.8871,5358.0820 4404,5473
2	4260.0000	5337.0000	1.7659598
3	4281.0000	5350.0000	2.3749254
4	4302.0000	5365.0000	3.1892015
5	4321.0000	5382.0000	2.6308863
6	4341.0000	5402.0000	1.7408107
7	4361.0000	5424.0000	0.56240958
8	4382.0000	5448.0000	0.23266472
9	4404.0000	5473.0000	0.0000000	Q4489.1994,5569.2638 4598,5651
10	4428.0000	5500.0000	0.45861062
11	4455.0000	5528.0000	0.55809286
12	4485.0000	5557.0000	0.53572793
13	4518.0000	5588.0000	1.3900942
14	4556.0000	5619.0000	0.62666202
15	4598.0000	5651.0000	0.0000000

始点 4237,5327 制御点 4306.8871 始終点 5358.0820 4404,5473 制御点 4489.1994,5569.2638 終点 4598,5651

山本ワールド

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

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

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

2次ベジェ曲線の近似がうまくいかないケースと対処法の計算例

うまくいかない近似例

うまくいかない近似例