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

icon 項目のみ表示/展開表示の切り替え
数式表示 \LaTeX(MathJax)

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


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

定数行列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になることも考慮すれば、楕円、放物線、双曲線、直線、二重線、点があるので、判別式で判定する。

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

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

座標リスト


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

係数行列K
1 2 3 4 5
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







制御点を表示

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

入力フォームへ移動

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

既知点からの距離
ixiyitbxbyd

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

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

うまくいかない近似例

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

既知点からの距離
ixiyitbxbyd
14237.00005327.00002.2204460e-16 4237.0000 5327.0000 1.2862197e-12
24260.00005337.00000.0049509460 4254.3045 5343.1157 8.3570896
34281.00005350.00000.010070127 4272.0346 5359.6272 13.155281
44302.00005365.00000.015529972 4290.7627 5377.0674 16.489332
54321.00005382.00000.021023723 4309.4176 5394.4386 16.996197
64341.00005402.00000.027182531 4330.1047 5413.7013 15.988366
74361.00005424.00000.033716199 4351.7897 5433.8921 13.516053
84382.00005448.00000.040803560 4375.0084 5455.5097 10.260529
94404.00005473.00000.048312213 4399.2620 5478.0894 6.9534855
104428.00005500.00000.056588981 4425.5852 5502.5941 3.5441459
114455.00005528.00000.065684698 4454.0152 5529.0580 1.4454296
124485.00005557.00000.075645895 4484.5525 5557.4808 0.65687467
134518.00005588.00000.086695866 4517.6962 5588.3264 0.44588397
144556.00005619.00000.098939356 4553.5216 5621.6636 3.6382684
154598.00005651.00001.0000000 4598.0000 5651.0000 0.0000000

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

うまくいかない近似例

上記の様なケースを検出するために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の要素の各項の値を算定
ixiyi誤差図形
14237.00005327.00000.0000000Q4306.8871,5358.0820 4404,5473
24260.00005337.00001.7659598
34281.00005350.00002.3749254
44302.00005365.00003.1892015
54321.00005382.00002.6308863
64341.00005402.00001.7408107
74361.00005424.00000.56240958
84382.00005448.00000.23266472
94404.00005473.00000.0000000Q4489.1994,5569.2638 4598,5651
104428.00005500.00000.45861062
114455.00005528.00000.55809286
124485.00005557.00000.53572793
134518.00005588.00001.3900942
144556.00005619.00000.62666202
154598.00005651.00000.0000000


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